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ABSTRACT 


Kacynski, Kenneth John. Ph.D., Purdue University, May 1994. Calculation of Propulsive 
Nozzle Flowfields in Multidiffusing Chemically Reacting Environments. Major Professor: 
Joe D. Hoffman, School of Mechanical Engineering. 

An advanced engineering model has been developed to aid in the analysis and design 
of hydrogen/oxygen chemical rocket engines. The complete multispecies, chemically 
reacting and multidiffusing Navier-Stokes equations are modelled, including the Soret 
thermal diffusion and the Dufour energy transfer terms. In addition to the spectrum of 
multispecies aspects developed, the model developed in this study is also conservative in 
axisymmetric flow for both inviscid and viscous flow environments and the boundary con- 
ditions employ a viscous, chemically reacting, reference plane characteristics method. 
Demonstration cases are presented for a 1030:1 area ratio nozzle, a 25 lbf film cooled noz- 
zle, and a transpiration cooled plug and spool rocket engine. The results indicate that the 
thrust coefficient predictions of the 1030:1 and the 25 lbf film cooled nozzle are within 0.2 
to 0.5 %, respectively, of experimental measurements when all of the chemical reaction 
and diffusion terms are considered. Further, the model’s predictions agree very well with 
the heat transfer measurements made in all of the nozzle test cases. 

The Soret thermal diffusion term is demonstrated to have a significant effect on the 
predicted mass fraction of hydrogen along the wall of the nozzle in both the laminar flow 
1030:1 nozzle and the turbulent flow plug and spool nozzle analysis cases performed. Fur- 
ther, the Soret term was shown to represent an important fraction of the diffusion fluxes 
occurring in a transpiration cooled rocket engine. 
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SECTION 1 
INTRODUCTION 

The computation of rocket nozzle performance is of critical importance to the rocket 
engine industry. Often, the rocket engine manufacturer needs to compare various engine 
operating systems before deciding on a specific configuration. Two items that are of criti- 
cal importance are nozzle performance and nozzle heat transfer characteristics. Often 
these two parameters are interrelated. This is most apparent in cases involving transpira- 
tion cooling or film cooling. Under these circumstances, the actively cooled boundary 
layer serves to thermally (and perhaps chemically) protect the nozzle wall from the 
exhaust gasses of the rocket engine core flowfield. In such conditions, the coolant allows 
the rocket engine to operate at conditions otherwise not attainable, such as high chamber 
pressures or low thrust levels. The engine designer must then ascertain whether the 
increased operating pressure or lower thrust level is worth the decreased performance of 
the engine due to active cooling. Consequently, a need exists for a comprehensive rocket 
nozzle performance prediction capability that can accurately determine nozzle perfor- 
mance and heat transfer characteristics for design and trade-off studies. 

A different need exists for individuals interested in improving rocket nozzle perfor- 
mance predictions and modelling. In these cases, it is important for the individual to ascer- 
tain the primary limitations of the nozzle performance prediction technique. A model 
which would include as many of the potential performance and heat transfer effects as 
possible and also allow the examination of individual effects, such as species diffusion or 
chemical kinetics, would be most desirable, with the additional constraint that the solution 
should be reasonably insensitive to the computational grid employed for the analysis. 
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1.1. Previous Efforts 

Rocket nozzle performance and heat transfer predictions have previously been 
accomplished primarily by the use of an inviscid core/boundary layer analysis technique, 
such as the JANNAF (Joint Army, Navy, NASA, Air Force) standard code, TDK/BLM 
(Reference 1), although others have used the full Navier-Stokes equations (References 2 
to 4). In either case, significant simplifying assumptions regarding the species diffusion 
and heat transfer mechanism are made. In general, all rocket engine performance predic- 
tion procedures known by the author have employed at least some of the following 
assumptions: 

• No chemical reactions 

• No thermal diffusion 

• No interdiffusional heat transfer 

• No diffusion thermo heat transfer 

• Laminar and turbulent Lewis numbers equal to unity 

• Simplified fluid transport coefficient models 

• Nonconservative evaluation of the finite difference equations 

• Extrapolated boundary conditions 

For many rocket engine applications, some or all of the above simplifications may be 
justifiable. In other cases, some, but not all, of the effects will be important. It is conceiv- 
able that there may exist conditions where all of the above processes should be included. 

1.2. Present Analysis 

In the present study, inclusion of the multispecies advection, production (i.e., chemi- 
cal kinetics), and diffusion terms, with full coupling in the total energy equation, was 
developed. Further, species and heat transfer by diffusion processes were modelled with 
consideration of the interrelation between species and heat transfer mechanisms (i.e., the 
Soret and Dufour effects). 

Modelling the multispecies rocket nozzle flowfield calculations was accomplished by 
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modifying the compressible, two-dimensional Navier-Stokes code, Proteus (Reference 5). 
In general, prior to these modifications, the Proteus code capabilities included the ability 
to numerically model the time-averaged Navier-Stokes equations for a single species per- 
fect gas (constant y). The Proteus code also included the ability to model a temperature 
dependent specific heat, thermal conductivity, and viscosity, fluid with the properties of 
air. Either an algebraic turbulence model or a two-equation (Chien k-e) turbulence model 
can be used for the closure of the Reynolds averaged turbulent transport terms. However 
the Proteus code can not adequately approximate turbulence conditions for a transpired 
wall and thus required modification. 

Axisymmetric flow can be modelled in the Proteus code although the analysis was 
nonconservative and could prove to be unstable (as described later). Modification of the 
code, based on a ‘conservation of one-dimensional flow’ concept was applied to both the 
inviscid and viscous terms to allow robust and accurate handling of the flow equations. 

Rocket nozzle performance prediction capabilities and accuracy were further 
enhanced by the inclusion of the following boundary conditions: 

• Characteristic equations at the boundaries 

• Mixed subsonic-supersonic exit conditions 

• Coupled temperature/species gradient wall conditions 

• Coupled heat transfer/transpired fluid wall conditions 

A very desirable attribute of the Proteus code is the available documentation (Refer- 
ences 5 to 7) and the modular program structure. As stated in the documentation “The Pro- 
teus two-dimensional Navier-Stokes computer program is a user-orientated and easily- 
modifiable flow analysis program for aerospace propulsion applications. Readability, 
modularity, and documentation were primary objectives during its development.” In main- 
taining the fundamental objective of the Proteus code, a new model, Tethys, was devel- 
oped for modelling chemically reacting and diffusing rocket nozzle flowfields. 
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SECTION 2 

GOVERNING EQUATIONS 

The governing equations for unsteady, viscous, compressible, multispecies laminar or 
turbulent flow in a two-dimensional axisymmetric coordinate system are quite compli- 
cated. They contain terms and correlations that are substantially less well characterized 
than the flow of a single species gas. While the general equations of multispecies fluid 
flow are reasonably well characterized for flows of engineering interest, the uncertainties 
of multispecies fluid flow lie in the determination of the chemical kinetic reaction rates 
and the transfer of heat, mass, and momentum by diffusion. 


2.1. Navier-Stokes Equations 

The general form of the axisymmetric equations governing the flow of a multiple spe- 
cies chemically reacting and diffusing gas can be written as follows (References 8 and 9): 


( 2 . 1 . 1 ) 


9(rQ) _ 8(rE) _ 9(rF) 
9t + 9x + 9r 


9(rE y ) d(rF y ) 
dx + 9r 


+ H 

v 


( 2 . 1 . 2 ) 

(2.1.3) 

(2.1.4) 


Q= [p, pu, pv, pe, pj] 


E = [pu, pu^ + P, puv, (p e + P )u, pju] 


F = [pv, puv, pv , (p e + P ) v, pjv] 


-,T 


,0, -ro. 


(2.1.5) 


r 9p 

H = 0, 0, r ^ 
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(2.1.6) 1 

T 

E = [0, x ,x , ux +vt — q_, — jj „ ] 

V xx xr xx xr 

(2.1.7) 

T 

F v = [°> VV“ t *r + Vt r,-‘lr.-ji,] 

(2.1.8) 

H v = [0,0,-X ee> 0,0] T 

For a chemical system comprised of n different species, the species continuity equa- 
tions (denoted by the i subscript) are solved for the density p. of species 1 to n-1, while p n 
is determined from the total density: 


n - I 

(2.1.9) Pn = P “ S P i 

i = 1 

2.2. Constitutive Equations 

The constitutive equations governing the diffusive transfer of momentum, energy, 
and mass, are somewhat complicated by the presence of multiple species fluid flow due to 
the coupled nature of the mass transfer and heat transfer terms. Mass transfer will be a 
function of the concentration gradients, temperature gradient, and the pressure gradient, 
and the species gradient, although for the particular applications of interest only the tem- 
perature and concentration gradients are of significance. Likewise, the heat transfer will 
also be coupled to the mass transfer, being comprised of conduction, interdiffusion, and 
Dufour heat transfer mechanisms, the latter two mechanisms being a function of the extent 
of mass transfer occurring. The general form of the constitutive equations is presented in 
Equations (2.2.1) through (2.2.8). 


( 2 . 2 . 1 ) 


Momentum transfer 


X 

XX 


_ du 
2p— + X 
ox 


du 1 d (rv) 

dx r 9r 


x 

TV 


dv 

2p.— — H X 
dr 


du 1 d (rv) 

dx r dr 


( 2 . 2 . 2 ) 
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(2.2.3) 

(2.2.4) 

(2.2.5) 

( 2 . 2 . 6 ) 


(2.2.7) 

( 2 . 2 . 8 ) 


x 

xr 



3u 3v 
3r + dx 


x ee = + x 


3u ld(rv) 
3x r 3r 


Heat transfer 


3T 

= " k §7 + X h . Ju 


n n 


MD 


i = 1 
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Diffusion thermo heat transfer 
(Dufour effect) 


Mass transfer 
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Thermal diffusion 
(Soret effect) 


j = V A ^ ^ — 

J ir Zj ij 3r T 3r 
j = 1 ' ' 

Ordinary Thermal diffusion 

diffusion (Soret effect) 
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The actual determination of the transport coefficients involved in momentum, heat, 
and mass transfer is very complicated, highly empirical, and rather inaccurate. This is 
especially true for the multicomponent species ordinary and thermal diffusion coefficients. 
The multicomponent solution for the coefficients of ordinary diffusion (Ay) and thermal 
diffusion (D?) are discussed in Appendix C. 


2 .3 , S pe cies S ourg g Fu n ction 

The species source function, a., is the term that represents the creation or destruction 
of a species by chemical reaction. Depending on the application, the term may be a rela- 
tively significant factor for both rocket engine performance (i.e., specific impulse) and 
heat transfer. 

For a mixture of n different gasses in which a total of m different chemical reactions 
can occur simultaneously, the general chemical reaction equation is discussed in Refer- 
ence 10 and is given below in Equation (2.3.1): 

n Kg n 

(2.3.1) XVi 

i = l K„j i=l 

where Kg and K b j refer to the forward and backward chemical reaction rates, respectively, 

/ n 

v and v refer to the stoichiometric coefficients of the reactants and products, respec- 
tively, and Aj denotes the specific chemical species. The equation governing the produc- 
tion or destruction of a chemical species is given by the following expression: 


(2.3.2) 


/ 




For the jth reaction. Kg is related to the backward reaction rate, K^j, and the equilibrium 
constant, K p j, by the following relation: 
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K f . _ -Av. 

(2.3.3) ^ = Kpj (RT) 1 

ft / 

where AVj = v. - Vj . 

For binary exchange chemical reactions, the backward reaction rate, K^, is a function 
of temperature only, and the functional relation is assumed to behave according to a modi- 
fied form of the Arrhenius equation (Reference 11): 


(2.3.4) 


a ; 


K bj = BjT J exp 


h_) 

rtJ 


For catalytic reactions (reactions involving a third body), a mean third body catalytic 
efficiency is used to determine a density weighted reaction rate: 


(2.3.5) 




y\ 

v 


where K bj - is a baseline backward chemical reaction rate for reaction j and the ratio 
Bjj/Bj is the catalytic efficiency of species i. Equations (2.3.4) and (2.3.5) readily show 
that the mechanisms for chemical reaction may be a strong function of species density. 
While the backward reaction rate of binary exchange reactions is a function only of tem- 
perature, catalytic reactions are an explicit function of the species density. Consequently, 
at low pressures, hence low species densities, chemical reactions tend to favor binary 
exchange mechanisms, whereas at high pressures (i.e„ combustor chamber pressure) cata- 
lytic reactions dominate. This is an important consideration in the design of many chemi- 
cal rocket engines. Since many of the chemical reaction mechanisms of interest involve 
catalytic reactions, rocket engine reaction efficiencies improve dramatically with increas- 
ing chamber pressures. A complete description of the species source term and its relation 
to the dependent variables ( p, pu, pv, pe, p. ) is given in Appendix D. 
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2.4. Thermal and Caloric Equations of State 

In general, the combustion gasses of a hydrogen-oxygen rocket engine are apprecia- 
bly in excess of the critical temperature of all of the individual species. Further, since the 
rocket engine gasses expand through the nozzle in a nearly isentropic fashion, the decay of 
fluid pressure maintains thermodynamic conditions that are significantly distant from the 
critical pressures and temperatures of all the chemical species. Consequently, the exhaust 
gasses of the rocket engine are assumed to obey the ideal gas law: 

( 2 - 4 -l) P = pRT 

The gas constant, R, is a weighted summation of the molecular weights of the individual 
species: 


(2.4.2) 


R 


= E S 


i = 1 


co. 

m: 


The assumption of ideal gas behavior also implies that the specific enthalpies and 
internal energies for individual species are solely functions of temperature and the average 
fluid enthalpy and internal energy are comprised of a mass fraction weighting of the indi- 
vidual species properties: 

T 

(2.4.3) = « iR + J C vi dT 

T R 

n 

(2.4.4) M = X ®i“i 

i = 1 

The fluid temperature can be determined by integration and rearrangement, and sub- 

2 2 

stitution of the definition of the total fluid energy (i.e., pe = p(«+ (u +v )/2)), 
which results in the following correlation: 
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pe - ^ (u 2 + v 2 ) - p £ (Oj (»iR~C vj T R ) 

(2.4.5) T = 

P X “Ai 

i = 1 

2.5. Time- Averaged Turbulent Flow Equations 

The presence of turbulence creates flow conditions that cannot be fully simulated 
with modem computational hardware. Consequently, the conventional technique of han- 
dling turbulent fluid flow is to time-average the fluid flow equations and model those 
terms that arise due to turbulent fluctuations. As the flow of concern may have appreciable 
Mach number variations throughout the turbulent region of the flow, a density-weighted, 
time-averaging technique has come into favor (Reference 12) as the preferred starting 


point for turbulence modelling. The complete derivation of the time-averaged fluid equa- 
tions and the assumptions associated with the averaging technique are presented in 
Appendix E. The resulting flow equations are presented below: 

3 (rQ) 

( 2 - 5 -D 3, 

d (rE) d (rF) _ 3 (rE v ) 3(rF v ) _ 

+ dx + dr +H dx + dr +Hv 

(2.5.2) 

T 

Q = [p, pu, pv, pe, pj] 

(2.5.3) 

T 

E = [pu, pu + P, puv, (pe + P)u,pju] 

(2.5.4) 

r _ __ -2 - _-|T 

F = [pv, puv, pv , (p e + P ) v, pjvJ 

(2.5.5) 

- r Bp - i t 

H — [o, 0, r 0, — 
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(2.5.6) 


E v = 


.T 


0, T + X 1 , X + X 1 , 
’ xx xx ’ xr xr ’ 


ux + vx -q +ux l + vx l -q* , - j; Y - i- 
xx xr n x xx xr n x ’ J i-x J i 


njt Jijc 


(2.5.7) 


Fv = 


-,T 


0, x + x l , x + x l , 

’ xr xr ’ rr rr ’ 


G \r + vx fr - q f + UX^ + vxf r - qj , - j i>r - j*. 


(2.5.8) 


H v [0,0, x Q0 x 00 , 0, 0] 


where the turbulent stress terms (denoted by the t superscript and defined in Appendix F) 
are assumed to have diffusive- like characteristics (i.e., the Boussinesq approximation) and 
can be defined as follows: 


(2.5.9) 

(2.5.10) 

(2.5.11) 

(2.5.12) 


x l = 

XX 


v £ + * 


3u 1 9 (rv) 
9x + r 9r 


x l =p l 

xr K 


Xlr = V ^ 


9u 9v 
9r + 9x 


rr 


x 00 - 2— v + X 1 


9u l9(rv) 
9x + r 9r 


9u ^ 1 9 (rv) 
9x + r 9r 


(2.5.13) 


9x 


- k ' I + I s . ji, 

i = 1 


(2.5.14) 


+ Xm!, 

i = 1 
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( 2 -5.15) jU=-^iM 

( 2 . 5 . 16 ) 

It is seen that the assumptions made concerning the transfer of momentum and con- 
duction heat transfer result in correlations that are entirely analogous to their laminar 
counterparts. As a result, the description of the turbulent fluid flow equations can be 
described in the same manner as laminar flow, with the exception that the constitutive 
equations for conduction heat transfer and for momentum transfer be represented by an 
effective thermal conductivity and viscosity terms, respectively. However, for multicom- 
ponent species transfer, no complete analogy exists between the laminar and turbulent 
energy and species transport terms, except in cases of binary species mass transfer with 
negligible thermal diffusion and negligible Dufour energy transfer. 

The introduction of turbulent viscosities, thermal conductivities, and diffusion coeffi- 
cients has been found to yield acceptable results for many engineering applications pro- 
vided the appropriate dependence of the turbulent properties on the fluid flow conditions 
and the geometry are adequately modelled. Typically, the modelling of the turbulent vis- 
cosity is performed by the use of either an algebraic or a differential equation turbulence 
model. In the Tethys model, only an algebraic (enhanced Baldwin-Lomax) turbulence 
model has been used. On the other hand, the modelling of the turbulent thermal conductiv- 
ity and mass diffusivities is performed by the introduction of turbulent Prandtl and Lewis 
numbers that are assumed to be either constant or an algebraic function of the turbulent 
viscosity and laminar Prandtl or Lewis numbers. A complete description of the models 
used in the Tethys fluid flow model and the enhancements required to model a multispe- 
cies rocket engine flow are described in Appendix F. 

Hereafter, the superscripts denoting time-averaged quantities will be discontinued 
(except in Appendices E and F) and it will be implied that the solution of the governing 
flow equations represent a time-averaged solution. 


5a>j 

dx 

3cbj 

dr 



13 


2.6. Nondimensionalization of Fluid Flow Equations 

In maintaining similarity with the original Proteus code, the multispecies Navier- 
Stokes equations are cast in a nondimensional form. The nondimen sionalizing factors 
used for the flow variables are listed in Table 1. 


Table 1. Nondimensional parameters used in the Tethys model. 


Nondimensional term 

Definition 

x°, r° 

x r 

L ref ' L ref 

u°,v° 

U V 

u ref u ref 

T°, p°, Pi ° 

T p Pi 

^ref Pref Pref 

P°, e°, R°, C p rcf 

P e p T ref r T ref 

2 ' 2 ' K 2 ' ^p 2 

Pref U ref Pref u ref u ref u ref 

t° 

u ref 

L ref 

V 

G. 

1 

a ref 

npO 

D ij 0 ' A u° D i 

D u , A 'J , 

^ab ref Pref ^ab ref Pref ^ab ref 

k°, \i° 

k eff Peff 
k ref Pref 
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The application of the above nondimensionalizing terms results in equations that are 
identical in form to the dimensional counterparts with the exception that various nondi- 
mensional groupings of reference terms are created. These groupings are listed in Table 2. 


Table 2. Nondimensional reference numbers. 


Nondimensional parameter 

Definition 

Reynolds Number, Re ref 

Pref u ref ^"ref 

Pref 

Prandtl Number, Pr fef 

Pref C p ref 

k ref 

Lewis Number, Le ref 

Pref^p ref ^ab ref 


k ref 

Source Number, S ref 

Pref U ref 

a ref L ref 


The Reynolds number (Re ref ) and Prandtl number (Pr fef ) are nondimensional 
groupings that are familiar to most engineers. The terms that arise due to multispecies 
fluid flow, namely the Lewis number (Le ref ) and the species Source number (S ref ) are 
less familiar terms. In fact, the Lewis number is occasionally defined as the inverse of the 
current definition (References 13 and 14) and, to the author’s knowledge, the species 
source number (or any term of a similar form) has never been used for chemical rocket 
engine studies. However, the source number is a very illuminating parameter. From this 
number, it can be seen that the species source term tends toward infinity when either the 
characteristic length (L ref ) or the characteristic species production/destruction term 
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( G re f ) tends toward zero. These two conditions are approached at the extreme ends of 
rocket engine propulsion systems. A small characteristic length usually represents a phys- 
ically small, low-thrust rocket engine. On the other hand, a small species production term 
may be a characteristic of a high pressure rocket engine which has essentially no dissoci- 
ated species present after combustion. For catalytic reactions, which are the reactions of 
greatest significance in a hydrogen-oxygen rocket engine, the species production rate 
increases with increased fluid density much faster than linearly, thus off-setting the refer- 
ence density term in the numerator of the species source term. Therefore, the species 
source term has a maximum significance somewhere between the two extreme engine 
sizes, the source number being large both at small engine sizes (frozen flow) and at large 
engine chamber pressures (undissociated flow). 

Table 3 illustrates the significance of the species source term for a spectrum of hydro- 
gen-oxygen rocket engines that are currently in production or development, including an 
evaluation of the Space Shuttle Main Engines operating at ten times the present opera- 
tional chamber pressure. In this table it is seen that the minimum source number occurs for 
engines in the size and pressure regime of the current Space Shuttle Main Engine, with 
satellite maneuvering systems being over an order of magnitude less sensitive to the 
effects of recombination chemistry. However, it is also worthy of note that in all cases 
examined, the source number does not exhibit the characteristic of being orders of magni- 
tude greater than unity, which would tend to lead one to conclude that chemical reactions 
should be considered in any chemical rocket engine analysis. However, the source number 
calculation is a reflection of the importance of chemical reactions to the species equations 
and is not a direct reflection of the importance of chemical reactions to engine perfor- 


mance. 
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Table 3. Species Source number for selected H 2 /O 2 rocket engines. 


Engine 

Thrust, lbf 

Chamber 

pressure, 

lbf/in. 2 

^ref 

10 X Space Shuttle Main Engine 

4.7 x 10 6 

3.3 x 10 4 

6.5 

Space Shuttle Main Engine 15 

4.7 x 10 5 

3.3 x 10 3 

3.5 

RL 10IIB 16 

2.0 x 10 4 

3.6 x 10 2 

7.3 

Space Station Maneuvering 17 

7.5 xlO 1 

7.5 xlO 1 

8.9 

Satellite Maneuvering 18 

5.0 x 10° 

2.5 x 10 1 

59.0 


The source numbers are based on approximate flow conditions at the throat using 
one-dimensional kinetics equations (Reference 19), using absolute value of all species 
production terms. 

Insertion of the nondimensional variables and reference terms into the multispecies 
Navier-Stokes equations and dropping the nondimensional superscripts and subscripts, the 
nondimensionalization of the Navier-Stokes equations yields the following equations: 


(2.6.1) 

3 (rQ) 3(rE) 3(rF) 8 ( r V a(rF v ) 

+ -. + H = -. + -. 

at ox or ox or 

(2.6.2) 

Q = [p, pu, pv, pe, pj] T 

(2.6.3) 

T 

E = [pu, pu 2 + P, puv, (p e + P ) u, PjU] 

(2.6.4) 

T 

F = [pv, puv, pv 2 , (pe + P)v, p i v] 

(2.6.5) 

a i T 

H = [o,0, r ,0, -r 


+ H. 
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( 2 . 6 . 6 ) 

(2.6.7) 

( 2 . 6 . 8 ) 


v ReL 


°' t xx’V UT xx + v ' c xr-pj>-^< Ji. x + j‘.x ) 


F = - 
v ReL 


°- V V u, xr + VX rr - K ' < ji,r + j'u ) 


H v = r e X tO. 0, - V 0, 0] 


The nondimensionalized constitutive equations have the form given by Equations 
(2.6.9) to (2.6.18) where it is recognized that the viscosity and thermal conductivity repre- 
sent contributions from both turbulent and laminar transport coefficients: 


(2.6.9) 

( 2 . 6 . 10 ) 


x = 

XX 


9u 


9u 1 9 (rv) 
9x + r 9r 


rr 


9v 

2ji— 4- X 
9r 


9u l9(rv) 
9x + r 9r 


( 2 . 6 . 11 ) 

( 2 . 6 . 12 ) 


x 

xr 


= ti 


9u 9v 
9r + 9x 


x ee = 2 »r +x 


9u l 9 (rv) 

9x r 9r 


(2.6.13) 



+Le X h i[ Ju+iL ] 

i = 1 


+ 


X 


MD 


co. 




Jj£ 

CO; 


J J^L 

CO 

J 
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(2.6.14) 


BT 


< lr = - k 3 j +U X h i[ Jij + Jij ] 


i = 1 


n n 

+ I x 


MD: 1 


CO. 


i = 1 j = 1 J i 4 U 


CO- < 0 - 

1 J 


(2.6.15) 


Jipc 


n 1 ^®j d [ aT 

^ A ij Bx T 9x 

j= 1 


(2.6.16) 

(2.6.17) 

(2.6.18) 




n-1 aoj. D T 9t 

A U dr T dr 
j= 1 


dco. 

;t _ _ (0 t l 

^iM a x 


dco. 

;t _ _ 0 t i 

4 iM dr 


The thermal and caloric equations of state have the same form as their dimensional 
counterparts. 



19 


SECTION 3 

GENERALIZED COORDINATE TRANSFORMATION 

The governing equations for axisymmetric two-dimensional flow presented in the 
previous section are not well suited for nozzle flowfield analysis due to the need for accu- 
rate boundary description and varying degrees of flowfield resolution desired throughout 
the flowfield. Consequently, coordinate transformations are often utilized in numerical 
computations to improve solution accuracy. This requires two major procedures; the trans- 
formation of the grid covering the physical domain of interest to a uniform orthogonal 
computational grid, and the transformation of the governing equations from the physical 
domain to the computational domain. The specific details of the coordinate transformation 
are contained in Appendix G, while the more significant details of the flowfield conserva- 
tion rules applied to the Tethys flow model are presented in this section. 

3.1. Generalized Transformation of the Advective and Pressure Flowfield Expressions 

Transformation of the inviscid portion of the governing equations results in the fol- 
lowing expression: 


(3.1.1) 


3(rQ) d(rE) 3(rF) . _ fviscousj 

3x 3r| v, terms , 

T 

[p, pu, pv, pe, p.] 

Q = 


(3.1.2) 


J 
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(3.1.3) 


E = 


£ x [pu, pu 2 + P, puv, (p e + P ) u, p.u] 

_ 

T 

£ r [pv, puv, pv 2 , (p e + P ) v, p.v] 


(3.1.4) 




Ti x [pu, pu 2 + P, puv, (pe + P) u, p.u] 


Ti r [j>v, puv, pv 2 , (pe + P) v, p.v] 


(3.1.5) H = 


0,0, r 


8(yVJ) 9(T] r P/J) 


- r c. 


3$ 


+ r 


dn 


, 0 , 


The metric derivatives and the inverse Jacobian ( J = 1 / J 1 ) can be evaluated from 
known geometric quantities: 


(3.1.6) 

(3.1.7) 

(3.1.8) 

(3.1.9) 


J 1 = Vr, - VC 

^X T| 

^ = — Jx 

r\ 

’I* = - Jr i 

ll r = JX 4 


(3.1.10) 
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3.2. Flow Conservation of the Advective and Pressure Terms 

A fundamental requirement for any numerical algorithm/grid transformation is to 
conserve flow properties. In a rectangular coordinate system, the flow conservation 
requirement dictates that the metric differencing applied to the flow variables be conserva- 
tive. An additional complication with axisymmetric flow, however, is that it is not entirely 
clear what flow properties are to be conserved. The reasoning and flow conservation 
requirements applied to the Tethys model is as follows: 

• One-dimensional flow properties in the axial direction (x) must be conser- 

vative 

• One-dimensional flow properties in the radial direction (r) must be conser- 

vative 

In general, the advective terms in the Navier-Stokes equations that contain deriva- 
tives in the axial direction have a form similar to the following: 

5(r<j)) 

(3-2.1) 

where, for example, <|> would be equal to pu in the continuity equation and pu + P in the 
axial momentum equation. Further, for one-dimensional flow. Equation (3.2.1) is equal to 
0. Consequently, integration of Equation (3.2.1) for the case of a one-dimensional fluid 
flow results in the following relationship: 

(3.2.2) r§ = C l +f(T) 

from which it follows that Cj must be equal to 0 and f(r) must be equal to n)> for the case 
where <|> is equal to a constant, which is the condition which exists in an axial, uniform 
flow environment. 

In a transformed coordinate system, Equation (3.2.1) can be divided by the Jacobian 
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and can then be written in the following strong conservation form: 


(3.2.3) 


i3(r«to = 

J dx 9r| V J J 


In order to satisfy the conservation of axial flow properties to a first-order require- 
ment (i.e., uniform flow), it is necessary that the following condition hold: 



A 



Physical space Computational space 


Figure 1. Discretization of flowfield about point 0 for an arbitrary grid. 


Consider now the arbitrary grid illustrated in Figure 1. Since central differencing of 
the flow variables is used in the Tethys code for both the flow variables and the grid met- 
rics, the following relation is required to maintain flow conservation in the axial direction: 


2Ar| 
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If central differencing of the grid metrics is applied and the radius terms are evaluated 
in an exact manner, the following discretized expression is obtained: 


r^<t>- 


(3.2.6) 


( r c" r b) 

2A£ 


~ r i<i>- 


( r d~ r a) 
2 At, 


( r c~ r d> 
2Ari 


-r 2 <b 


< r b~ r a ) 

2Aq 


2At\ 2A t, 

Axial flowfield formulation using exact radius technique 


To be conservative in the uniform, axial flow case. Equation (3.2.6) must be equal to 
zero when <J> is equal to a constant. Apparently, the above discretization does not satisfy 
the requirement of flowfield conservation (by being equal to 0) for an arbitrary grid in the 
case where the dependent flow variable, <J), is constant. In the Tethys code, flow conserva- 
tion is satisfied by using a second-order accurate approximation for the radius terms at 
Nodes 1 to 4: 


(3.2.7) 


r l = 


r d + r a 


t 2 = 


r K + r 
b a 


r 3 = 


r b + r c 


r 4 = 


r c +r d 


Substitution of the above relation into Equation (3.2.5) yields the following expression: 


(3.2.8) 


(r d- r a>. 

4A4 * 4A4 * 


2Ati 2AE, 

Axial flowfield formulation using average radius technique 


4Ari ^ 4At| 


It is seen that the use of radius averaging results in axial flow conservation when the 
dependent variable is constant for all arbitrary grids. 

In the radial conservation scheme, the variation of the independent variable r must be 
considered. Unlike the axial derivative terms where the original presence of the radius 
inside the derivative was misleading to the one-dimensional flow case, the metric term is 
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crucial in the radial flow terms. Therefore, the first order conservation of radial flow stipu- 
lates that the following condition must be true: 


(3.2.9) 


3 




* *r* 


= 0 when r<j> is constant 


, J ) dev J 
General radial flow conservation requirement 


For example, the parameter <j> is equal to pv in the continuity equation. Additionally, the 
maintenance of a constant pressure flowfield is crucial and requires that 


(3.2.10) 



+ 



0 


when 



For both the radial flow and the pressure terms, no unconventional handling of the 
metric terms is required to maintain flowfield conservation. The exact radial formulation 
outlined in Equation (3.2.6) is adequate for Equation (3.2.9) and Equation (3.2.10). How- 
ever, it must be noted that the above procedure of flowfield conservation is unable to 
exactly satisfy conditions existing in the radial momentum equation due to the weak con- 
servation form of the radial pressure gradient. 


3.3. Generalized Transformation of Diffusive Terms in the Navier-Stokes Equations 

The transformation of the diffusive terms in the Navier-Stokes equations follows the 
same trend as the advective terms with the additional requirement that the diffusive fluxes 
of momentum, heat, and mass transfer must also be transformed. A complete description 
of the transformation is contained in Appendix G. The resultant expressions are presented 
below. 


(3.3.1) 


Advective and pressure terms = 


d(rE v ) 

~ar 


3 (rF v ) 

an 


+ 



where the transformed viscous terms for the continuity, x-momentum, r-momentum. 



25 


energy, and species equations are 


r 

(3 ' 3 ' 2) ^'r^jL 0, u \x + v ' C xr-‘lx /Pr . ~ Le <j«+j!,x ) /Pr ] 

+ RtT7 |_°’ ’ l ’ T ’ Ut xr + v V-<ir /Pr - - Le O.j+Jij ) /Pr 

(3.3.3) F v = Aj [0, t xx , t xr , ut xx + vt xr - q x /Pr, -Le (j u + j}, ) /Pr" 

r ... 1 T 

+ Sri L 0 ’ X xr ’ Z " ’ “ T *r + vt rr - % /Pr - Oij + Ji, ) / p r 


(3.3.4) 


H v = p 7 


1 [ 0 , 0 > -t ee ,°, 0 

J L _ 


The fluxes can be expressed in terms of a transformed coordinate system as illustrated 
below, where the viscosities and thermal conductivity represent both laminar and turbulent 
contributions. 


( 3 - 3 .5) t xx - 2*i £ xgr + Tl 


„ du 3u 1 . r„ du du 1 ( Y d(rv) d(rv)V 

+ T ^xdq + ^xd^ + 1 ^x5rj + 7 di; +T *r’dri ^ 


(3- 3 -6) !„ = 


. f dv dv 1 f du du 1 / d(rv) d (rv) 

2p ^55- + p r5^r +x ^xg^ +r| xg^ + 7nt - 55 — + T 'r‘g^ 

J L 


r „ du du „ dv dv 

(33 ' 7) X xr = ^ ^ + fl x ^ 
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(3.3.8) 


, r e 3T dl “I _ r , , 

q r = “ k W + ’'r 3^ + Le 2 h i [ Ji/ + h 


3T 


as 


dT\ J 


i = 1 




n n 


+ LeR ^ I I 


MD 


CO- 


tiTuTi J L 


Jj£ _ 

CO. CO. 
I J 


(3.3.9) q = 


■ - k [ S* 5? + ’•x 3^ ] + Le I h i [ iu + iL ] 

i = 1 


n n 


+ LeR n 2 2 


MD 


p 2* 2* 

K i= lj = 1 J i 4 iJ L 


Ji,x _ Jj^c 

CO. CO. 

i J 


(3.3.10) hj 


(3.3.11) Ju 


(3.3.12) 


(3.3.13) 


n - 1 


j=l 


\ 


3co. 8co. I Dj r 9 t 8T “I 

^ as + ^ aif J ~ t|_ as +T1 r _ 


n- 1 


j=l 


A 'j 


3co. 3 co. 

£ + 0 -5-^ 

3S '* 3rj 


D 


t 3T 3T 

3S + dr\ 


Jij “ ^iM 


3c0j 3c0j 

lf + T1 r J 


Ji^t ^iM 


dcOj 3cOj 

S£ + n x ^ 


In transformed coordinates, the nonconservative shear stress term, x 06 , can be written as 


(3.3.14) x ee = 2 [L- + 1 


„ 3u 3u 

^ 3f + ^ 3ri + 


a (rv) 'N 


1 ( 3 (rv) 3 (r\ 

7^-35- + 


The relationships between the metric derivatives S x > T| x , £ p T| p the Jacobian, J, and 
physical space derivatives are presented in Section 3.1. 
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3.4. Conservation of Diffusion Terms 

Similar to the conservation requirements for one-dimensional inviscid flow, it is 
highly desirable to develop schemes that can retain one-dimensional diffusion solutions 
independent of the grid transformation. The methodology in which this was accomplished 
follows. To the author’s knowledge this formulation represents the first conservation con- 
cept applied to viscous terms. 

To exemplify the procedure applied to one-dimensional viscous flow conservation, 
consider the flow of a fluid contained between a rod moving at some constant velocity 
within a stationary pipe, as illustrated in Figure 2. Due to its analogy to planar flow, the 
flow condition will be termed axisymmetric Couette flow. 





Figure 2. Row between concentric tubes. 


In this flow situation, the shear stress and velocity distribution are governed by the 
following expression: 


(3.4.1) 


d(rc ) 

v xr 7 

3r 


= 0 


where the shear stress, x xr , simplifies to the following: 
(3.4.2) T xr = ^ 


The solution of this flow situation is easily obtained by integration and substitution. The 
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shear stress has the following radial distribution: 


(3.4.3) 



The velocity distribution is obtained in a similar manner. 


(3.4.4) u = C 1 ln(r)+C 2 

where the viscosity has been assumed to be constant. 

Consider now the possibility of obtaining the exact flow solution numerically and 
with an arbitrary grid. This requires that the following expression (obtained by transfor- 
mation of Equations (3.4.1) and (3.4.2) and dividing by the Jacobian) be satisfied when the 
velocity profile given by Equation (3.4.4) is used: 



The evaluation of the transformed, viscous term outlined in Equation (3.4.5) can be 
performed about a point in an arbitrary grid, as illustrated in Figure 3. 



Physical space Computational space 


Figure 3. Discretization about point 0 for an arbitrary grid. 

A central difference discretization of Equation (3.4.5) can be obtained by averaging about 
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midpoints, as shown in Equation (3.4.6): 



Consider now the substitution of the exact value of the viscous shear stress term, 
for the case of axisymmetric Couette flow, that is, Equation (3.4.3) and apply central dif- 
ferencing to the metric derivative terms. The resulting expression that needs to be satisfied 
for numerical conservation of the shear stress term is: 


x c~ x d x 3 -x l x 3 _x 

4A^Ari C 4A^Aq C + 4A^Aq'“ " r 4A^Ar| 


1 X h~ X a 

ic + . C 


x c” x b 


x 4“ x 2 


+ _c + - 

4ArjA^ 4AqA^ 


x 4 - x 2 X. x 

c - - -A-^c = o 


4AqA^ 


4At|A^ 


By careful inspection, it is readily observed that the above formulation is satisfied. 

The remaining requirement is that the exact solution of the discretized shear stress 
terms can be realized for an arbitrary grid. This requires that the following condition be 
true: 

(3.4.8) V = 7 = (*5 = n| r ^[C,ln(r)+C 2 ] + Hl) r JL[C,ln (r) + C 2 ] 


As an example, consider the evaluation of the shear stress term at node 4 in Figure 3. 
Using first-order differencing of the velocity gradient in the £, direction and central differ- 
encing in the ti direction, the following requirement for flow conservation are obtained: 
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(3.4.9) 


C _ p C ln ( r 4 /r o) , „ C ln ( r c /r d) 
r 4 ~ + 2At| 


Inspection of the above expression reveals that it is necessary that the evaluation of the 
metric terms contain terms containing logarithms. Substitution of the metric derivatives in 
terms of physical space derivatives and rearrangement yields the following requirement 
for shear stress conservation: 


(3.4.10) = 

r 4 


ln(r 4 /r Q ) 


+ x 


l n ( r c /r d ) 

5 2At| 


A solution which satisfies the above conservation requirement involves the evaluation of 
the metric terms and inverse Jacobian in the following manner: 


(3.4.11) 


x n = 


(*c~ x d> ( r 4~ r 0> 
2At| r 4 ln(r 4 /r Q ) 


(3.4.12) 


x t = 


(*4 ~ V < r c ~ r d> 
2A^r 4 ln(r c /r d ) 


(3.4.13) J = x 5 r n - 


( x c" x d> ( r 4"V 
2A^Arj 


The metric derivatives are comprised of a standard backward or central difference opera- 
tor and then ‘tuned’ by a correction factor. The correction factor is approximately equal to 
1 + Ar/(2r) , therefore retaining the accuracy of the differencing technique and being 
significant only in regions near the centerline or in regions of extremely coarse grid reso- 
lution. 

Substitution of the expressions for the inverse Jacobian and spatial difference metrics 
into the flux conservation requirements of Equation (3.4.10) result in the following 
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expression for the case when the shear stress is inversely proportional to the radius: 

(X 4 -X 0 ) (r c -r d ) (X c -X d ) (r 4 -r 0 ) 

<3 ' 4 ' 14) r 4 2A^Ar| r 4 2 A^At| 

< x c ~ x d> < r 4 “ r 0> < x 4 - x <p ( f c " r d> 

r 4 2A^Ari + r 4 2A^All 

It is readily observed that the above equation is satisfied and the flux term is therefore fully 
conservative in the one-dimensional flow situation. 

It is also necessary to insure that the flux terms that should be zero in a Couette Flow 
situation are indeed equal to zero. For instance, the shear stress X xx should be zero. How- 
ever, only careful handling of the metric terms and the introduction of correction factors as 
performed on the terms insures that this is indeed the case. 

The technique applied to the metrics that act on scalar quantities, like temperature, is 
somewhat less obvious, as it would be desirable to have the following conditions of flux 
transport satisfied: 


(3.4.15) 


3t 

dn _ 

5x “ 


CjA 



whenT = Cjln(r) + C 2 x + C 3 


Unfortunately, the exact satisfaction of the above equation is not entirely clear. However, 
applying the same methodology to the scalar quantities as applied to the velocities, the fol- 
lowing conservation criteria are readily obtained: 


3t 

Tt 


= Cj/r 



Ar 


when T = C j In ( r) + C 2 x + C 3 and — — > 0 


(3.4.16) 
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For regions near the centerline, the condition of axial flux conservation is relaxed and 
the conservation of radial fluxes is enforced: 

(3.4.17) ^ = Cj/r when T = Cjln (r) +C 3 

In general, the rules that are applied in the Tethys code for viscous flow conservation 
satisfy the following criteria: 


• The Differencing of the metrics is never outside of the computational box. 

The box outlined in Figure 3 contains all of the necessary locations for 
evaluation of the dependent and independent variables. 

• The thin shear layer assumption cannot conserve diffusive fluxes. Except in 

cases where the grid transformation is one-dimensional in nature, the 
application of the thin shear layer assumption removes any possibility of 
diffusive flow conservation. This is because the balancing of the shear 
stress and velocity profiles requires the use of all grid metrics and differ- 
ence equations in order to obtain an expression that is independent of the 
grid transformation applied. Therefore the thin shear layer approximation 
is not recommended. 

• Radial correction factors are applied to all terms except the radial velocity. 

Since the primary diffusive transfer of radial momentum will be in the 
axial direction, conservation of axial momentum transfer dictates that no 
correction factor be applied. 
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SECTION 4 

BOUNDARY CONDITIONS 

In the original Proteus code, boundary conditions were ill-adapted for most nozzle 
flow situations. At the nozzle inlet, second-order extrapolation of velocities was the rec- 
ommended technique for solving the velocity profile. Even worse, at the nozzle exit, com- 
plete extrapolation of all the fluid properties was recommended for supersonic flow 
whereas the recommended technique for subsonic flow was to specify the pressure and 
extrapolate the velocities and temperature. More discouraging was the fact that the Proteus 
code did not have the means of distinguishing between the two conditions and required 
that the user specify that the flow at any exit grid point is either subsonic or supersonic. 
Quite obviously, in the case of a mixed subsonic/supersonic exhaust as experienced in real 
rocket engines and in transient simulations, the exhaust boundary condition was totally 
unacceptable. 

In the present study, the Tethys code used a reference plane characteristics method. 
This technique was found to provide a good representation of the actual flow physics and, 
for inflow and outflow boundaries, it was well adapted for implementation into the struc- 
ture of the original Proteus code. 

Wall boundary conditions in the original Proteus code required modification for the 
added complexity of multispecies fluid flow, with the additional need for including a cou- 
pled fluid injection/heat transfer boundary condition. 

In the sections that follow, a discussion of the enhanced boundary condition capabili- 
ties developed in the present study for the evaluation of hydrogen/oxygen rocket nozzles 
are presented. 
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4.1. Constant- Ti Reference Plane Method of Characteristics 

Appropriate specification of inflow and outflow boundary conditions has been shown 
to have a significant effect on the overall accuracy and convergence rate of the flowfield 
equations (Reference 20). The method of characteristics has long been used as the pre- 
ferred procedure for determining the appropriate equations that need to be solved at fluid 
flow boundaries (References 21 to 23) due to its ability to decompose the propagation of 
information into information carried by wave surfaces (the envelope of which is a Mach 
conoid) and information carried along a pathline, as illustrated in Figure 4. In the constant 
reference plane characteristic technique applied in this study, the assumption is made a 
priori that the characteristic lines of interest will lie solely in one direction and the applica- 
ble compatibility equations in the specified direction are developed. In the case of a com- 
putational scheme employing transformed coordinates as used in this study, the most 
obvious approach for obtaining characteristic and compatibility relations is to derive 
expressions valid along a constant-eta (t|) or constant-xi (^) plane, as illustrated in Figure 
5. 
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Figure 4. Characteristic surfaces and lines for unsteady two-dimensional flow. 


x 



Figure 5. Constant-T] reference plane. 
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The resulting relations for the reference plane characteristic and compatibility equa- 
tions for an inviscid, nonreacting fluid are given below. A detailed derivation which 
includes the viscous and chemical reaction terms applied in the Tethys model is presented 
in Appendix H. The results for the constant- T| reference plane are presented in this sec- 
tion. The resulting expressions for a constant-^ reference plane are analogous (requiring 
only the switching of the £ and Tj coordinates) and are also applied in the Tethys code. 


Table 4. Constant- rj reference plane method of characteristics. 


Characteristic Equation 

Compatibility Equation 

^ = U 
dx 

dP - a dp = \y 4 dx 

dx 

M u- ^x dv = <5rV2-W dT 

«.u 

dx 

dco i = ^4 + i dT 


+ ^dP + pa^du + pa^dv 


= [ ( a 2 V J + v 4 ) J^ + ^ + pa (S x V 2 + ^ r V 3 ) ] dx 

mi 

fel + ^ dP_ P a ^ x du_ P a ^r dv 


= [ (a 2 Vi + V 4 ) Jt* + ^-pa (^ x v 2 + ^ r V 3 ) ] dx 


In Table 4, the contravarient velocity in the ^-direction, U, is defined as U = u£,^ + v^. 
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In these results, the inviscid source and cross-flow terms, \|/j to \y 4 + . , for chemi- 
cally frozen, inviscid fluid flow are given below: 

P v du dv dp 

V 1 " T " ^ ^ ~ V ^ 


¥ 2 - 


= - V 


du 


dP 


v 3 - 


v h. 

drj 


\ dP 
P dri 


»4--v5 + .*V* 


^4 + i 


= - V 


d<a. 

Si* 


where the contravarient velocity in the r| -direction, V, is defined asV = ur| x + vr| r . 

The characteristic and compatibility equations are of the same form in either dimen- 
sional or nondimensional terms. 

The characteristic equations are extremely powerful in that they determine the appro- 
priate flow equations that need to be solved at computational boundaries. In regions where 
the characteristics lie outside of the computational domain, the corresponding compatibil- 
ity equations must be replaced with appropriate boundary conditions. For instance, at a 
subsonic nozzle inlet with inflow, as illustrated in Figure 6, the only characteristic that lies 
within the computational domain is the d^/dx = U-a ^ characteristic. The 
compatibility equations applicable on the characteristic curves which are outside the com- 
putational domain (shown in the figure as the shaded region) must be replaced by appro- 
priate boundary conditions. Similar requirements exist at the fluid outflow boundaries and 
along the walls. 
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4.2. Subsonic Nozzle Inlet Boundary Conditions 

As discussed in the previous section, only one characteristic curve lies within the 
computational domain of a subsonic inflow nozzle. Consequently, boundary conditions 
must be specified to form a complete set of inlet conditions. In a chemical rocket engine 
analysis, one normally has knowledge of the total specific enthalpy of the gasses, the mass 
concentration of each species (usually determined from an equilibrium chemistry analy- 
sis), the flow angle, and the stagnation pressure. An alternative formulation that is equally 
valid for both chemically reacting and nonreacting flow is to specify the fluid entropy 
instead of the stagnation pressure. 

4.3. Nozzle Exit Boundary Conditions 

The exhaust from a propulsive nozzle is primarily supersonic at steady-state condi- 
tions. In these circumstances, all of the fluid characteristics lie within the computation 
domain and no additional boundary conditions need to be specified. However, during tran- 
sients and at regions very near the wall, the flow may be subsonic. Further, many verifica- 
tion experiments and analyses are performed at subsonic conditions, thus requiring the use 
of a subsonic nozzle exit boundary condition. For subsonic nozzle exhausts, one charac- 
teristic lies outside of the computational domain, and the standard boundary condition to 
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impose at the exit is the specification of the fluid static pressure. During transients, of 
course, the ability to discern between subsonic and supersonic flow is essential and is a 
feature of the present analysis. 

4.4, Nozzle Centerline and Free Slip Wall Boundary Conditions 

At a nozzle centerline and on an inviscid free slip wall, it is highly desirable to model 
those conditions where a gradient between the nozzle boundary and the adjacent fluid flow 
point in the interior of the computational flow exists. This can be achieved by applying 
one waveline characteristic equation and the following boundary conditions: 

• The normal velocity, V n , equals zero. 

• The entropy gradient, 3 s/dn, equals zero. 

• The stagnation enthalpy gradient, 3h 0 /3n, equals zero. 

• All species gradients, 3co./9n , equal zero. 

In light of the fact that the above boundary conditions (i.e., the last three listed) would 
be better modelled with additional characteristic relations that apply along the fluid path- 
line (which unfortunately were ill-suited for the structure of the numerical procedure 
employed, as discussed in Appendix A and Section 5), the additional approximation is 
made on the stagnation enthalpy, entropy, and species concentrations, that the gradient in 
the physical normal direction is equivalent to the gradient in the computational normal 
coordinate system. 


4.5. No Slip Wall Boundary Conditions 

At a no slip wall, the boundary condition dictates that both the radial and axial com- 
ponents of velocity be zero, 
u = 0 
v = 0 

Further, there can be no flux of any species normal to the wall. Therefore, the following 
species flux condition, assuming that no pressure driven species fluxes are of concern, 
must be satisfied: 



40 


(4.5.1) 



n * d(0 i 

X A ij^r 

j = i 


DfdT _ ]i l Le l d(B i _ 

T dn p r t chi ^ 


Note that only in cases of negligible thermal diffusion will the gradient of the mass frac- 
tions be equal to zero. 

Another boundary condition that needs to be specified at a no slip solid wall is the 
thermal boundary condition. A general linear case is specified by the mixed boundary con- 
dition: 


(4-5.2) T w = aq w + b 

For an isothermal wall boundary condition, a = 0, while for a specified heat flux boundary 
condition, b = 0. For a solid wall boundary, since there can be no species flux, the heat flux 
equation reduces to the conduction equation: 

(4.5.3) q = - k 

H w a n 

The remaining boundary condition that needs to be specified at a no slip solid wall is 
usually obtained by assuming that the normal pressure gradient is equal to zero. This is 
usually an excellent steady state approximation, even in cases of transpired flow. How- 
ever, instead of specifying a boundary condition that is essentially valid only at steady 
state, a more robust boundary condition can be obtained by employing one waveline char- 
acteristic equation. 


4.6. Transpired Wall Boundary Conditions 

Transpiration cooling adds considerable complexity to the appropriate boundary con- 
ditions that are necessary for accurate description of the conditions existing at the wall. 
However, a few boundary conditions are similar to a solid wall. First, in order to be con- 
sidered a solid boundary, the tangential velocity along a transpired wall must be zero. 
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(4.6.1) V t = 0 

Similar to the no slip solid wall case, and somewhat more applicable, one characteristic 
waveline equation is applied at the wall. 

The remaining boundary conditions that need to be applied to the transpired wall are 
less apparent and somewhat dependent on the application of interest. Figure 7 represents a 
conventional transpiration cooled nozzle wall. 


%’ (P v >„- h w 


! 


: i I Z J 1 

Coolant flow, ( P v) Q , 

Figure 7. Transpiration cooled wall. 

A control volume analysis of the transpired wall yields the following relation: 
(4.6.2) (P v ) 0 h o = [<P v >„ h + <ln] w 



where the subscript w indicates conditions existing at the boundary between the transpired 
wall and the fluid flowfield and the subscript o indicates conditions existing upstream of 
the transpired wall. Since mass diffusion may occur at the boundary of a transpired wall, 
the heat flux will be composed of conduction, diffusion, and Dufour energy transfer com- 


ponents: 

(4.6.3) 
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The species density boundary conditions are also complicated somewhat by the tran- 
spired wall boundary condition since some species diffusion can occur and the appropriate 
boundary condition is to relate the known species flux at the wall to the diffusion and 
advective transport terms: 

(4.6.4) h(pv) n + jj, n + it, ] w = [ pveflj ] o 

Two flow conditions of interest may exist. In the ‘analysis mode’, knowledge of the 
normal flow rate from the transpired region is known. In these cases, Equations (4.6.2) and 
(4.6.4) provide the necessary boundary conditions for determining the wall species fluxes, 
the heat fluxes, and the wall temperature. On the other hand, a ‘design mode’ condition is 
often more desirable. In this case, the wall temperature is specified and Equations (4.6.2) 
and (4.6.4) are used to determine the total mass flow rate necessary to maintain the wall at 
the required temperature. 

The boundary conditions presented for the transpiration cooled wall are in dimen- 
sional form. Nondimensional forms will contain reference Reynolds number, Prandtl 
number, and Lewis number expressions. 
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SECTION 5 

OVERALL NUMERICAL ALGORITHM 

The numerical algorithm used to solve the governing fluid flow equations in the 
Tethys model employs a scheme developed by Beam and Warming (Reference 24). In this 
scheme, the implicit finite difference operators (see Appendix A) are divided into opera- 
tors that act in the £, direction, those that act in the T| direction, and those that act on the 
transient term. In the procedure used by the Tethys model, the implicit operators represent 
central-difference operators analogous to the operators that act on the explicit side of the 
equation. This is in contrast to other numerical methods that employ linearized flux-split- 
ting methods on the implicit side of the equation and central difference techniques on the 
explicit side of the equation (References 25 and 26). 

The interior flowfield model is capable of marching forward in time by means of 
either Euler implicit or explicit methods, or by marching forward in time using a trapezoi- 
dal implicit or a three-point backward implicit method. The Euler methods are first-order 
accurate in time while the other procedures are second-order accurate in time. The Euler 
implicit procedure was selected as the most appropriate methodology to use throughout 
the course of the demonstration and verification studies because only the steady state solu- 
tion was of interest. 

The solution methodology used by the Beam and Warming scheme involves the pro- 
cess of first solving the flowfield equations that lie on a constant-q line by using an 
implicit operator that acts on the transient and the fully coupled £ direction terms. This 
procedure is then repeated along all of the interior constant-q lines using a highly vector- 
ized solution methodology. The procedure is then repeated along the interior constant-^ 
lines with the explicit side represented by the changes in the flow properties that were 
observed in the q direction sweep. In the Proteus model, the constant-^ boundary points 
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were then updated (after an rj and a £, line sweep) to maintain a more rigorously correct 
transient solution. However, as the cases of interest in these studies all involved steady 
state solutions, the boundary points were all updated at the same time after the entire inte- 
rior solution had been marched forward one time step. 

Since only the interior computational lines arc swept through in an implicit manner, 
boundary conditions cannot be applied implicitly except for those terms that lie on the 
same interior computational line corresponding to the sweep direction. Likewise, since the 
boundary points at the comers of the flowfield are not located on any interior computa- 
tional line, they need to be updated in a purely explicit manner. 

The CFL (i.e., Courant, Friedrichs, and Lewy, Reference 27) number was used as the 
basis for the maximum time step that could be taken as the solution proceeded from an ini- 
tial profile towards a steady state condition. In the Tethys model, the CFL number is based 
purely on the time that it takes a Mach line to travel from one grid point to an adjoining 
point. This is in contrast to other computational models that require a viscous correction to 
the CFL number to account for the propagation speeds of the viscous terms (Reference 
25). Because the Tethys code used an implicit formulation of the diffusion terms, a viscous 
correction was not deemed to be necessary. 

In the majority of the verification and demonstration cases studied, the flowfield solu- 
tion was obtained by marching every grid point at a local CFL number that was in the 
range of 0.7. However, in some cases, somewhat faster convergence was obtained by 
using a mixed CFL number criteria, developed by Knight (Reference 28), where the max- 
imum local CFL number in the q direction was set equal to a time step value of around 4.0 
and the maximum local CFL number in the % direction was not allowed to exceed a value 
in the region of 0.4. In yet other cases, a global CFL number criteria was used to examine 
the accuracy of the implicit formulation. In these cases the CFL number varied from max- 
imum values of 10 to as high as 10000, depending on the application. 

In general, the initial flowfield used to initiate fluid flow problem should employ a 
property distribution that is as close to the steady state solution as could be obtained. In 
many cases this can be accomplished by using one-dimensional flow solutions as a start- 
ing point. However, this was not always practical and some cases were started by using a 
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semisteady truncated nozzle solution. The initial profile for these truncated geometries 
varied from a one-dimensional chemical equilibrium solution to a near zero velocity con- 
dition. Once a reasonably close solution for this geometry was obtained, the remainder of 
the nozzle was added on to complete the analysis. 

In general, a nozzle flowfield problem that started from an initially low velocity dis- 
tribution would prove to be unstable as the boundary conditions applied at the exit of the 
nozzle would impose a static pressure that was too great a change from the initially stag- 
nant properties. Significant alleviation of this difficulty was obtained by scaling the maxi- 
mum change in pressure at the exit as being equal to 15 % of the product of the static 
pressure multiplied by the local CFL number. 

In nearly all of the demonstration cases examined, the heat transfer rate was the last 
parameter of engineering interest to reach an acceptable level of steady state convergence. 
Consequently, convergence of the nozzle flowfield was determined to occur when the tran- 
sient variation in heat flux to the nozzle wall was deemed to be negligible. The flowfield 
predictions from the Tethys model were inspected at periodic intervals to ascertain when a 
converged solution had been obtained. In general, the computational resources required of 
the problem and the accuracy to which a heat transfer prediction was desired was used as 
the guide as to when an acceptable convergence of the heat transfer rate was established. 
In all cases, the minimum convergence criteria occurred when the absolute value of the 
local heat fluxes at selected points had less than a 2 % change after 1000 iterations. In 
many other cases, the convergence criteria was significantly more stringent. 

In general, a flowfield problem required on the order of 10,000 iterations to achieve 
the desired level of convergence. The computational resources required for this level of 
convergence varied from 2 to 50 hours of computational time on a Cray Y-MP computer. 
The chemically frozen, binary species flowfield simulations required the minimum com- 
putational effort in the viscous analysis cases performed, while the maximum computa- 
tional resources were required when the Tethys model was operating at its maximum 
flowfield simulation capacity of six species and eight chemical reactions. 
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SECTION 6 

VALIDATION STUDIES 

Extensive validation studies were performed to verify the prediction accuracy of the 
Tethys code since there were significant changes made to the Proteus code to obtain the 
Tethys code. Validation of the following effects was deemed critical: 

• Conservation of one-dimensional flow 

• Reference plane characteristics boundary conditions 

• Chemical kinetics 

• Transpired fluid boundary conditions 

• Multicomponent diffusion (ordinary and thermal) 

• Multicomponent heat transfer (conduction, interdiffusion, and Dufour) 

• Turbulence model validation 

• Fluid thermodynamic and transport properties validation 

With the exception of the first two validation requirements presented above, the vali- 
dations are either contained in the Appendix or verified within a reasonable level of cer- 
tainty in the demonstration studies or in other unreported validation tests. Further, the 
solution procedure of the Tethys fluid flow model is fully implicit in nature with complete 
coupling between the species and fluid flow equations, including a complete description of 
the dependency of the fluid pressure, temperature, and total energy on species concentra- 
tion. A discussion of the solution methodology used throughout the validation studies is 
given in Appendices A and B. 

The two validation studies that follow (conservation of one-dimensional flow and ref- 
erence plane characteristics boundary conditions) are given more extensive coverage since 
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the application of these concepts is both a unique and an interesting aspect of the Tethys 
code. 


6.1. Conservation of One-Dimensional Inviscid Row 

Reformulation of the manner in which the original Proteus code analyzed the radius 
in the axisymmetric flow equations was required. As discussed in Section 3, the original 
Proteus code evaluated the radius exactly, but one-dimensional axial flow conservation 
shows that the radius in the axial flow terms are conservative when the radius is evaluated 
in terms of averaged quantities. About the rj point: 


( 6 . 1 . 1 ) 


_ r ^ + i,ti + 1 +r ^+ i,n — i 
VM 2 


( 6 . 1 . 2 ) 


_ r £+l,T|+l +I ^-1,T| + 1 
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and likewise for other points. 

The ‘nonconservative nature’ of using the exact formulation for the radius can be rea- 
soned to be a second-order effect since the difference between the exact and the averaged 
radius procedure is also a second-order difference. Therefore, for cases where the second- 
order derivative of the radius with respect to the grid metrics is zero, no error in using the 
exact evaluation of the radius occurs, and the flow is fully conservative. Consequently, 
verification of the accuracy of a finite difference procedure using metrics that have sec- 
ond-order derivatives that are zero is misleading in the more general case. 

The nonconservative error associated with formulating the pressure using the exact 
radius and a source term formulation is of the following nature: 


1 9(rP) p 3(r£ r P/J) 3(rri r P/J) 

J dr J + dr| ^ 


Substitution of the expressions for the metrics and inverse Jacobian yields the following 
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condition for pressure conservation: 

d(rx P) 3(rxtP) 

(61 ' 4) + — 5ii P (x ^ " r $V = 0 

when the pressure, P, is equal to a constant 

If central differencing of the grid metrics and fluid pressure is used, the conservation 
requirement given in Equation (6.1.4) is satisfied only for cases in which any one of the 
axial metric derivatives, such as x^, is equal to zero. This condition is common for many 
rocket nozzle types of analysis, although requiring that this metric be equal to zero for 
pressure conservation is ill-advised since some nozzles flow cases are much more effi- 
ciently executed in an environment in which all of the grid metrics are nonzero in nature. 
Likewise, for cases in which the radial derivative is a function of only one variable, which 
is typical for many pipe flow simulations, no solution error exists in the above formulation 
of the pressure. 

Verification of the nonconservative nature of the original metric terms and the conser- 
vative nature of the reformulated terms will be shown for inviscid flow in a pipe with two 
different grid generation cases. In one case the metrics are a function of only one indepen- 
dent variable, in which case both the exact and the radial-averaging technique should yield 
stable solutions. In the other case, the metrics are a modest function of both x and r and 
consequently the exact radial formulation would be expected to yield nonconservative 
solutions that may not converge. 


6.2. Inviscid Pipe Flow - One-Dimensional Grid Transformation 


The 11X11 grid used to demonstrate the nature of the flow in a simply packed grid is 
shown in Figure 8 where the grid packing involved simple one-dimensional power rela- 
tionships: 


( 6 . 2 . 1 ) 
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Figure 8. One-dimensional transformed grid for inviscid pipe flow calculation. 

Characteristic boundary conditions, as described in Section 4, were used at all bound- 
aries. The back pressure was set to 14.69 lbf/in. 2 , while the stagnation pressure at the inlet 
was set to yield a Mach number of 0.4. The fluid was air (75.83 % N 2 , 23.13 % O 2 , and 
1.03 % At by weight) and the stagnation enthalpy of the mixture was zero. The exact solu- 
tion (to 5 significant digits) was used as the starting point for all of the test cases exam- 
ined. No artificial diffusion was applied. The local CFL number was set to 0.7. 

Three discretized flowfield approximations were considered. One formulation, identi- 
fied in Figures 9 and 12 as the standard technique, employed an exact radius evaluation 
and formulated the pressure terms in the radial momentum equation as a combination of a 
strong conservation term and a source term, as shown in Equation (6.1.3). The second 
flowfield approximation method, identified in Figures 9 and 12 as the exact radius tech- 
nique, formulated the pressure in a weak conservation form, while evaluating the radius in 
an exact manner. The third finite difference approximation method used the complete con- 
servation concepts outlined in Section 3, and is identified in Figures 9 and 11 as the con- 
servative technique. 
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The results for the inviscid pipe flow simulation for all three flowfield formulations is 
presented in Figure 9. It is evident from Figure 9 that all of the evaluation techniques con- 
verge. The standard formulation converged the slowest, while the exact radius and conser- 
vative techniques have almost identical convergence rates. 



Figure 9. Convergence rate of standard, exact radius, and conservative techniques for 
inviscid pipe flow with a one-dimensional grid transformation. 


6.3. Inviscid Pipe Row - Two-Dimensional Grid Transformation 


A pertubation of the previous grid description was made for the case of the general- 
ized grid transformation. In this case, the grid metrics were determined using the follow- 
ing relations: 


(6.3.1) 
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Figure 10 presents the resulting grid generated when the above relations are applied. 
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Figure 10. Two-dimensional grid used in the inviscid flow conservation study. 


The exact same boundary conditions, initial conditions, and time stepping procedure 
was applied to the evaluation of inviscid flow in a pipe with a two-dimensional grid trans- 
formation as was applied to the test cases performed with inviscid pipe with a one-dimen- 
sional grid transformation. These boundary conditions include a back pressure of 14.69 
lbf/in. 2 , a stagnation pressure at the inlet set to yield a Mach number of 0.4, and an inlet 
stagnation enthalpy of zero. The fluid was air and the local CFL number was set to 0.7. No 
artificial diffusion was applied. 

Three test cases were examined. The cases examined are described in Section 6.2 and 
represent various finite difference approximations to the radius and the radial pressure gra- 
dient. The three different test cases are identified as the standard technique, the exact 
radius technique, and the conservative technique. 

In Figure 11, the stable nature of the conservative technique is apparent as the con- 
vergence of the solution continues at about the same rate as that observed for the one- 
dimensional transformed grid case, with the maximum change in density continuously 
decreasing as the number of iterations increases. 
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Figure 1 1 . Convergence rate of conservative technique for inviscid pipe flow with a 
two-dimensional grid transformation. 

The standard and the exact radius techniques proved to be unstable when the grid was 
changed to a more general two-dimensional case, as illustrated in Figure 12. The noncon- 
servative nature of the formulations are apparent as the standard formulation becomes 
unstable in less than 250 iterations, while the exact radius formulation begins to diverge 
after about 1000 iterations and becomes unstable in less than 1500 iterations. 



Figure 12. Divergence rate of standard and exact radius techniques for inviscid 
pipe flow with a two-dimensional grid transformation. 
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6.4. Couette Flow 

The potential benefits of applying a flow conservative formulation to the diffusive 
terms was examined for the cases of flow between two parallel plates and flow between 
two concentric cylinders. The 6X11 computational grid used for both studies is two- 
dimensional and nonorthogonal to the wall as presented in Figure 13. The bottom surface 
was set to a velocity of 1000 ft/sec and the upper surface was stationary. The exact veloc- 
ity profile was imposed at the inlet. The fluid examined was air and the Reynolds number 
based on the bottom surface velocity and gap between the fluid surfaces was equal to 0.1. 



Figure 13. 6X11 Computational grid employed for the analysis of Couette flow. 
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The results for the planar Couette flow study are given in Figures 14 and 15. On the 
overall scale shown in Figure 14, it appears as if both the conventional diffusion equation 
technique and the conservative formulation used in the Tethys model are very accurate. 
However, closer examination of the actual velocity error as presented in Figure 15 indi- 
cates the superiority of the conservative technique. It is observed that the velocity error 
exceeds 5 percent in the conventional analysis while the flow conservative analysis yields 
an order of magnitude less error in the velocity. 
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Figure 14. Predicted velocity profiles for a Couette flow simulation. 

While it might be expected that the Tethys formulation would yield exact results for 
this flow situation, this is not the case as the modelling of the inviscid term pu 2 was 
designed to be conservative for the case where the axial momentum was independent of 
the y-location. In the case of Couette flow, the term is not constant throughout the flow- 
field. Consequently, the introduction of a nonzero gradient of axial momentum in the 
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radial direction resulted in a nonzero contribution to this advective term in the axial 
momentum equation. 



u = 1000 ft/sec q Grid point # u = 0 ft/sec 

Figure 15. Axial velocity error at the computational midline of a Couette flow 
simulation. 


The results for the annular axisymmetric Couette flow simulation (which used the 
same grid as used in the planar flow case but was displaced such that the ratio of the inner 
radius to outer radius was equal to ten) are given in Figures 16 and 17. The trends in the 
results are similar to the planar flow condition with both the standard finite difference for- 
mulation and the conservative Tethys formulation having significantly more error than 
their planar flow counterparts. The maximum error observed with the conventional, or 
standard, formulation exceeds 10 percent, which is nearly an order of magnitude greater 
than the maximum error observed with the Tethys viscous conservative scheme. Again, 
the failure of the Tethys modelling scheme to exactly yield the correct solution is attrib- 
uted to the nonconservative nature of the discretized advection term in the axial momen- 
tum equation. 
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Figure 16. Predicted velocity distribution for axisymmetric Couette flow simulation. 



u = 1000 ft/sec r\ Grid point # u = 0 ft/sec 

Figure 17. Axial velocity error at the computational midline of an axisymmetric 
Couette flow simulation. 

It is worthy of mention that in all of the Tethys flow cases examined above, the con- 
servation of species mass fraction for all of the constituents of air (nitrogen, oxygen, and 
argon) was completely satisfied. 
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6.5. Supersonic Radial Source Flow 

Radial source flow represents a one-dimensional flow condition that needs to be con- 
sidered prior to making judgement of the two-dimensional accuracy of a numerical flow 
simulation. For the purposes of this study, it represents a flow situation that allows the 
examination of how well the Tethys model conserves species concentrations, stagnation 
pressure, stagnation enthalpy, and mass. Further, it provides a baseline case from which 
one can compare chemically reacting solutions, as is shown in Appendix J. 

Figure 18 pictorially presents a general case of radial source flow. Solutions for this 
type of flow are usually not determined at the centerline but instead the results are scaled 
from some arbitrary starting radius as represented by the dashes in Figure 18 with the 
understanding that some type of source flow could produce those conditions at the starting 
line. 
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Figure 18. One-dimensional radial source flow. 


The inlet boundary conditions applied to the source flow simulation are presented in Table 
5, which corresponds to an inlet Mach number of 1.5 and species mass fraction corre- 
sponding to a stoichiometric mixture ratio of hydrogen and oxygen. Characteristic bound- 
ary conditions were employed at the exit while reflective boundary conditions were used 
along the upper and lower walls. 
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Table 5. Inlet boundary conditions used for frozen source flow examination. 


Property 

Inlet Value 

Temperature 

6060.6 R 

Density 

0.005 lbm/ft^ 

Velocity 

7171.9 ft/sec 

“h 2 

0.0156 

“° 2 

0.0829 

“h 2 o 

0.7500 

“oh 

0.1050 

“h 

0.0025 

“o 

0.0190 


The 4X21 grid used to demonstrate the accuracy of the radial flow terms in the 
present metric and flow methodology is shown in Figure 19. The grid represents a radially 
packed grid that should provide a good check of the ability of the code to model one- 
dimensional radial flow. 

As can be seen in Figure 19, the radial grid packing was one-dimensional and of the 
form outlined in the following expression: 

r = r inlet + [ r Exit~ r inlet ] 


£-1 


^Max 1 
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Radius, ft 

Figure 19. Computational grid used for supersonic planar source flow simulation. 


The nonlinear second-order artificial diffusion model (see Appendix N) was applied. 
The magnitude of the diffusion coefficient was examined at values of 0.5 and 0.125, which 
is twice as large and one-half as small, respectively, as the recommended value of 0.25 
(Reference 29). Since no attempt was made to alter the artificial diffusion model to force 
one-dimensional flow conservation, its effect on the conservative nature of the flow vari- 
ables required examination. 

In Figure 20, the predicted and exact Mach number distribution for the flow situation 
described above are shown. The discrepancy between the exact results and those predicted 
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by the Tethys flow model are hardly distinguishable, which in a rather global sense dem- 
onstrates the modelling accuracy of the radial derivative terms contained in the Tethys 
flow model. However, it can be observed that the artificial diffusion has a distinct effect on 
the numerical predictions, over predicting the Mach number at low diffusion values and 
under predicting the exit Mach number at higher levels of diffusion. 



Figure 20. Mach number predictions for a chemically frozen 100:1 area ratio planar 
source flow calculation. 

The extent to which the mass flow rate, stagnation enthalpy, and stagnation pressure 
are conserved is presented in Figures 21 to 23, respectively, for the two different levels of 
diffusion examined. It is readily observable from the mass flow rate and stagnation 
enthalpy conservation cases that artificial diffusion has a deleterious effect on the conser- 
vation of mass and stagnation enthalpy. On the other hand, the artificial diffusion has no 
clear effect on the stagnation pressure and in fact the higher diffusion level actually gave 
better conservation of stagnation pressure than the lower diffusion level. All of these 
trends might be expected since the evaluation of the transformation metrics were designed 
such that one-dimensional radial flow would be exactly conservative in mass and energy 
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with the exclusion of the consideration of any possible nonconservative effects of the arti- 
ficial diffusion terms. Quite obviously, complete conservation of mass flow and stagnation 
enthalpy is not realized and the lack of conservation of these variables is clearly a direct 
function of the level of artificial diffusion applied. On the other hand, conservation of stag- 
nation pressure requires that all of the flow equations be completely satisfied. However, 
the modelling of the radial pressure gradient term was performed such that it was equal to 
zero for the case where the pressure was constant throughout the flowfield. For one- 
dimensional radial flow, the pressure can vary significantly throughout the flowfield and 
can consequently result in a lack of conservation of stagnation pressure, in addition to the 
nonconservative contributions of the artificial diffusion. 
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Figure 21. Mass flow rate conservation study for a 100: 1 planar source flow simulation. 



Figure 23. Stagnation pressure conservation study for a 100:1 planar source flow simulation. 
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6.6. Inlet Exit, and Wall Reference Plane Characteristics Verification 

Verification of the constant-^ reference plane characteristics technique for nozzle 
inlets and exhausts is demonstrated throughout most of the verification and demonstration 
tests performed. In this particular example, special emphasis is placed on the fundamen- 
tals of the subsonic inlet, the supersonic exit, the inviscid wall, and the centerline tech- 
niques, leaving more advanced verification tests that involve chemical kinetics, diffusive 
transport, and artificial diffusion, to the more difficult demonstration tests. 

The implementation of the constant-T| reference plane characteristic solution at any 
time level is performed as shown in the Figure 24. The flow variables and appropriate 
equations are determined by linear scaling with conditions that would exist at the maxi- 
mum explicit time step, dx M0C . For example, the solution of the compatibility equation, 
dP - a dp = \j/ 4 dx, along the characteristic curve d^/dx = U, is performed in the fol- 
lowing manner (using the nomenclature shown in Figure 24): 


(6.6.1) dP - a 2 dp = \|/ 4 dx - P?+/ - P? 



- v 4 Ax MO c 


where the terms without subscript and superscript nomenclature are evaluated at time 
level n and the average of the values existing at spacial levels i and i+1. Values at time 
level n+1 can be decomposed into quantities at time levels n and a delta form, as shown in 
Equation (6.6.2): 


( 6 . 6 . 2 ) 


pn + 1 
r i + 1 
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i + 1 


+ 


d? n 
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P n 

r i + 1 


+ AP n 


+ 1 


The expression for the density has a similar form. 

Substitution and rearrangement of Equation (6.6.1), using the concepts shown in 
Equation (6.6.2), result in an discretized expression valid along the constant-T| plane char- 
acteristic pathline. 
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(6.6.3) AP" + , -a 2 Ap" + j = " 


P n P n 

r i+l r i 


+ a “ 


P n , - p" 
K i+1 


For time steps less than the time step corresponding to the characteristic curve, the 
solution is scaled by employing a similar triangle methodology, as illustrated in Figure 24. 
The resulting general compatibility equation (for time steps less than or equal to the time 
step along the characteristic curve) is given in Equation (6.6.4): 


r... 


(6.6.4) AP? +1 -a^A Pi " +1 =[ ¥4 -J? _ (Pf +1 -Pf) 


+ a ‘ 




While this technique is only first-order accurate in time, a desirable attribute to this 
method is that the steady state solution of the compatibility equations are independent of 
the size of the time step used in obtaining a solution since the interior flow point algorithm 
also contains only a global Ax operator on the explicit side of the equation (see Appendix 
A), similar to Equation (6.6.4). 
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Figure 24. Numerical implementation of the constant-q reference 
plane method of characteristics. 
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Verification of the characteristic technique for rocket nozzle flowfield analysis will be 
performed using the JPL nozzle tested by Cuffel et al. (Reference 30). The nozzle geome- 
try is shown in Figure 25, while the 45X11 grid used for the flowfield analysis is presented 
in Figure 26. This grid scheme is reasonably close to the grid applied by Driscoll (Refer- 
ence 20) whose results will be used for comparison between the reference plane character- 
istic method and the bicharacteristic technique. The boundary conditions applied for this 
nozzle analysis case are presented in Table 6. 
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Figure 25. Schematic of small radius of curvature nozzle tested at JPL. 
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Figure 26. 45X11 grid used for JPL transonic nozzle flow study. 


Boundary 

Boundary conditions 

Inlet 

Constant-r) reference plane 
P D = 70 lbf/in 2 
T 0 = 540 R 
v = 0 
Dry air 

Exit 

Constant-Tj reference plane 
9co. 

55 1 = 0 

Wall 

Constant-^ reference plane 

v n = o 

£-0 

0T1 

aP o 

3-° = 0 

9ti 

9g). 

5— 1 = 0 

Centerline 

Constant-^ reference plane 
v = 0 

£•- 

8co. 

3- 1 = 0 

9ri 
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In Figure 27, a comparison of predicted Mach number contours, using the reference 
plane characteristic method, the bicharacteristic method, and experimental measurements 
is shown. The numerical predictions of the reference plane characteristic method are rep- 
resented by dashed lines, while the solid lines represent predictions using the bicharacter- 
istic method. When followed along contours originating at the centerline, the shaded 
contours, solid lines, and experimental measurements, are all near the same origination 
point from which one can track subsequent disparities in the interior flowfield. 

Excellent agreement between experimental measurements and numerical predictions 
using the present reference plane characteristic method is observed, especially in the tran- 
sonic region. In the supersonic region, the prediction of the compression wave that origi- 
nates at the tangency point between the downstream circular arc and the conical nozzle is 
modelled reasonably well along the wall, although the flow is extensively diffused in the 
interior of this region. Further, it is observed that the slight streamline curvature between 
the nozzle centerline and the adjacent interior grid points is predicted in the supersonic 
region. In general, the results of the reference plane characteristic method appear to pre- 
dict results at least as well as the bicharacteristic method. 



Bicharacteristic method 

Reference plane characteristics method 



Figure 27. Comparison of Mach number contour predictions and experimental measurements for the JPL nozzle. 
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SECTION 7 

DEMONSTRATION STUDIES 

The cumulation of a diverse set of concepts and models intended to ultimately 
improve the predictive capabilities for a chemical rocket engine is demonstrated in this 
section. For illustrative and demonstrative purposes, three unique hydrogen-oxygen 
chemical rocket systems are examined. First, a high area ratio (1030:1) rocket nozzle will 
be analyzed and discussed in detail. The future application of such a nozzle is envisioned 
to include orbital transfer, lunar, and planetary missions. The uniqueness of the high area 
ratio nozzle is that it represents a very high performance rocket nozzle, and it has an 
extensive experimental database that includes performance, heat flux, and static wall 
pressure measurements. Further, the nozzle is one with which the author is very familiar, 
having been involved in the reporting of the experimental data. A second hydrogen-oxy- 
gen rocket engine of current interest is a low thrust (25 lbf) engine that is being envi- 
sioned for satellite thruster applications and for the orbital maintenance of a man-tended 
space station. The uniqueness of this nozzle concept is that it is extensively film cooled 
with hydrogen, and it operates at an overall stoichiometric mixture ratio while having an 
oxygen/hydrogen mixture ratio of 20 in the core of the engine. This nozzle concept has 
also been extensively tested. Engine performance, global heat transfer, and static wall 
pressure measurements are available for comparison. The final configuration studied is a 
transpiration cooled plug and spool rocket engine. The uniqueness of this configuration is 
that a discrete section of the rocket engine has hydrogen coolant transpired along the 
wall. Hydrogen coolant mass flow rate and wall temperatures were measured in this 
experiment 

A complete discussion of the nozzle geometries, operating conditions, modelling 
techniques and assumptions employed, and a comparison of experimental results and 
numerical predictions follows. 
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7.1. 1030:1 High Area Ratio Nozzle 

The high area ratio hydrogen/oxygen rocket nozzle illustrated in Figure 28 (Refer- 
ences 31 to 33) was tested at the NASA Lewis Research Center. Gaseous hydrogen and 
oxygen were injected into a combustion chamber that had a subsonic area ratio of 4.22. In 
the chamber, the propellants were burned to near completion, with the combustion effi- 
ciencies in each test ranging from 93 to 99 percent The flow then expanded through the 
converging-diverging nozzle and was exhausted into a diffuser where the pressure was 
partially recovered. Following the cooling of the exhaust gasses by a water spray, the non- 
condensable gasses were entrained by nitrogen driven ejectors and exhausted into the 
atmosphere. The diffuser inlet and the complete rocket nozzle assembly were contained 
within a vacuum test capsule, as illustrated in Figure 29. 
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Figure 28. High area ratio nozzle test program conducted at NASA LeRC. 
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Figure 29. Rocket engine test facility used for evaluation of 1030: 1 area ratio nozzle. 


Experimental measurements of the nozzle thrust and mass flow rate were made. Fur- 
ther, wall static pressure and temperature measurements were made in the supersonic 
region of the nozzle. Heat fluxes were computed by employing the transient response of 
the nozzle thermocouples. 

A detailed uncertainty analysis of the experimental program (Reference 34) indicated 
that the instrumentation uncertainties involved in the nozzle performance measurements 
were less than 1.3 %. However , the potential biasing effect of throat contraction or 
expansion during the testing of the nozzle was not accounted for. 

The validity of assuming that axial conduction and radiation heat transfer were negli- 
gible in the measured heat transfer rates was examined during the experimental program 
(Reference 32) and estimated to amount to less than a 2 % error. An additional 
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investigation of the experimental data (performed for this study) indicates that the nonlin- 
ear variation in the temperature response of the thermocouples and the use of locally con- 
stant material properties in the data reduction may result in an additional data uncertainty/ 
bias of 3 %. 


7.1.1. Analytical Assumptions 

For analytical comparison, several important assumptions were made concerning the 
environment under which the rocket nozzle performed. The major assumptions regarding 
the flowfield included complete mixing at the nozzle inlet, the magnitude of the exhaust 
pressure at the nozzle exit, and laminar flow throughout the nozzle. 

The chemical composition at the inlet of the nozzle section was assumed to be fully 
mixed and in chemical equilibrium, leaving the lack of complete mixing as a ‘correction 
factor’ that will be applied when engine performance is compared to experimental mea- 
surements. Furthermore, the existence of the trace species of monatomic and diatomic 
oxygen were neglected since the flow was very fuel rich and predictive results from a one- 
dimensional equilibrium flow analysis indicated that the effect of these trace species was 
negligible, resulting in less than a 15 degree R change in the stagnation temperature of the 
nozzle and less than a 0.1 % difference in specific impulse. 

The exhaust pressure at the nozzle exit is an important parameter as the nozzle could 
experience some level of ‘boundary layer feedback’ where the subsonic flow region of the 
boundary layer is affected by the ambient pressure at which the nozzle operates. In turn, 
this can affect the pressure on the nozzle wall for some distance upstream of the nozzle 
exit plane. Since the nearest static pressure tap on the nozzle was located 2.7 inches from 
the exit, a small extrapolation (4.4 %), using the interior nozzle static pressures was 
required to establish the back pressure at which the nozzle was operating. It should be 
noted that this back pressure was 30 % lower than the capsule pressure, thus indicating 
that the nozzle was somewhat overexpanded. 

The nozzle flowfield was also assumed to be laminar throughout. Extensive studies 
and correlation with experimental and analytical predictions indicate that the rapid accel- 
eration of the exhaust gasses in the throat region of the nozzle Taminarizes’ the flowfield 
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(see Reference 32). Additional studies have indicated that the laminarized profile is stable 
throughout the supersonic region of the nozzle (Reference 35). Since the transonic and 
supersonic regions of the nozzle are the area of greatest concern, the assumption of lami- 
nar flow throughout the nozzle is reasonably justified. 

7.1.2. Boundary Conditions 

The boundary conditions applied in the Tethys code are presented in Table 7. The ref- 
erence plane characteristics method (MOC) was employed along all of the boundaries, 
including the no-slip wall boundary. The wall temperature distribution employed for this 
analysis, illustrated in Figure 30, was based on experimental measurements and interpola- 
tion when possible, and was based on engineering judgement in regions outside of the 
experimentally measured area. The discontinuous temperature slopes observed near the 
throat region exist due to the application of a ceramic coating in this region only. 


Table 7. Boundary conditions applied to the LeRC 1030:1 nozzle. 


Inlet 

Exit 

Centerline 

Wall* 

P 0 = 360.0 lbf/in. 2 

MOC* 

MOC 

MOC 

h 0 = -15 Btu/lbm 

MOC 

3h 0 /3r| = 0 

u = 0 

v = 0 

MOC 


v = 0 

MOC 

MOC 

5s/3t| = 0 

T = T (x) 

Chemical equilibrium* 

MOC 

dcOj/dri = 0 

ji = 0 

* The inlet mass 
fractions are: 
go H2 = 0.093302 
cd h = 0.003052 
a) 0H = 0.020253 

“h.o = 0883393 

*At subsonic 
contravarient 
velocities, static 
pressure is 
specified: Pgxit = 
0.0264 lbf/in. 2 


* For a free slip 
wall, the 
boundary 
conditions used 
are discussed in 
Section 4.4. 
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Figure 30. Wall temperature distribution for the high area ratio rocket nozzle analysis. 


7.1.3. Computational Grid Used for Analytical Study 
In Figures 31 and 32, the 91 X 51 computational grid used for the analysis of all three 
viscous test cases is presented. The inviscid flow analysis used a slightly more coarse 
computational grid, as discussed in the following paragraphs, but was nearly indiscernible 
to visual inspection. 



Every third rj and t, line printed 



Axial length, in. 


Figure 31. Computational grid employed for high area ratio nozzle analysis. 



Radius, in. 



Figure 32. Throat region view of computational grid employed for high area ratio nozzle analysis. 
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The most noteworthy aspect of the computational grid is that the centerline length is 
significantly shorter than the nozzle length. This ‘centerline coring’ of the nozzle takes 
advantage of the fact that the supersonic flow can only propagate information within a dis- 
crete region (the domain of influence) and test runs with coarse grids indicated that the 
selected grid includes the complete domain of influence. The net result of this ‘coring’ in 
the centerline region was that it allowed significantly more grid points to be placed in 
physical locations where they would be of use to resolve the flowfield and the conditions 
that exist along the wall. 

As illustrated in Figure 32, significant packing was applied along the nozzle wall and 
in the throat region. The grid is also slightly tighter packed along the nozzle centerline, 
and the grid packing along the wall is much tighter in the subsonic/transonic regions than 
near the nozzle exit. 

The grid was generated by first determining the location of the £, end points that 
would lie along the nozzle and along the centerline. This was done by the application of 
the transformation discussed in Anderson, Tannehill, and Pletcher (Reference 36), Equa- 
tions (5.221) and (5-223), which are repeated below: 


(7.1.1) 


x = x 


Throat 


1 + 


sinh [t(^-B)] 
sinh (tB) 


where B is defined as follows: 


(7.1.2) 


B = 


T" ln 

2t 


1 + 


1 + 


— T 


f X Throa| N 

) V X Total ) 

^Throat 
X Total j 



The packing factor, T, was set equal to 8.0 along the nozzle wall and was equal to 6.0 
along the nozzle centerline. For the centerline transformation, the total centerline length 
was equal to twenty percent of the actual nozzle length. 

The transformation of the T) grid lines was performed according to Equation (5-216) 
in Anderson, Tannehill and Pletcher (Reference 36). 
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(7.1.3) 




1 

^Total 



— *n 



l-Tl 

+ 1 


where 1 is the distance along a computational line to the centerline, and the coefficients a 
and b are related to the packing factor, p, by the following relations: 

(7-1-4) a = p+ 1, b = P - 1 

For the viscous analyses cases (i.e.. Test cases 2 to 4), the packing factor, p, was defined as 
a function of ^ through the following relation: 

(7.1.5) p = 1.001+ (1.01 - 1.001) £? 

For the inviscid analysis, the packing factor, p, was set equal to 1.01. 

The grid points established using the above grid transformation procedure were over- 
written in the region near the centerline (the nearest seven points) to allow somewhat 
greater packing near the centerline. The algebraic transformation applied to every grid 
point number, N, in this region is given below: 


(7.1.6) 



N - 1 
8-1 


1.5 


-°- 5 Urr 


Here <j> is equal to the fractional distance from the centerline to the eighth point located 
away from the wall and along a constant ^ line. 

7.1.4. Analytical Study and Results 

Four different Tethys model simulations were performed to examine the significance 
of the modelling assumptions that were employed. These different analysis methodologies 
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are presented in Table 8. Here it is seen that the analysis comprised a spectrum from chem- 
ically frozen, inviscid flow to chemically reacting diffusing flow with thermal mass trans- 
fer and Dufour effects included. 


Table 8. Various analysis methodologies employed for 1030:1 LeRC nozzle. 


Tethys 

analysis 

Test case 1 

Test case 2 

Test case 3 

Test case 4 

Chemical 

kinetics 

No 

No 

Yes 

Yes 

Viscous 

flow 

No 

Yes 

Yes 

Yes 

Ordinary 

diffusion 

No 

No 

Yes 

Yes 

Thermal 

diffusion 

No 

No 

No 

Yes 

Interdiffusion 
heat transfer 

No 

No 

Yes 

Yes 

Dufour 
heat transfer 

No 

No 

No 

Yes 


7 . 1 .4. 1 . Mach Number Contour Predictions 

In Figure 33 to Figure 40, the predicted Mach number contours for all four test cases 
are shown. The comparison between the inviscid and viscous cases (Test cases 1 and 2, 
respectively) in Figures 33 to 36 indicates that the Mach number at the exit for the inviscid 
simulation is slightly higher than the Mach number in the viscous analysis, which is an 
indication of the fact that the boundary layer, because of its displacement thickness, will 
tend to reduce the effective area ratio of the expansion, with a resultant decrease in the exit 
Mach number. Comparison between the frozen and kinetics throat region cases (Test cases 
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2 and 3, respectively) in Figures 35 to 38 reveals that the Mach number in the chemically 
frozen analysis is somewhat higher throughout the supersonic portion of the nozzle than in 
the chemically reacting analysis. Physically, this is as expected since the chemical reac- 
tions release sensible energy and tend to slow down the rate of expansion, a ‘Rayleigh 
line’ effect. Finally, the last two test cases which are shown in Figures 37 to 40, which are 
distinguished only by thermal diffusion and Dufour heat transfer effects, have extremely 
similar Mach number profiles. 

Close inspection of the inviscid analysis case shown in Figures 33 and 34 indicates 
that the flow near the wall of the nozzle in the supersonic region of the flow undergoes a 
large amount of compression. The result of this compression in the supersonic region of 
the nozzle is that the Mach number contours are nearly parallel to the wall. 

The Mach number contours along the centerline in the supersonic region of the flow 
are retarded. This is possibly because of the computational grid employed for the analysis, 
as the constant £, lines go from left -to- right as the physical location goes from the center- 
line to the wall. On the other hand, the actual wavelines striking the center of the nozzle 
emanate from the leftward direction. As a result, the reference plane characteristic method 
uses an adjacent point that is well outside of its domain of dependence, and the predictive 
accuracy of the centerline boundary point would be expected to suffer. Fortunately, the 
analysis of the centerline region appears excellent in the transonic and slightly supersonic 
region, and the error that occurs near the end of the centerline region has little effect on the 
wall solution. 
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Figure 37. Mach number contour predictions for the LeRC high area ratio nozzle using a viscous, chemically reacting 
analysis (Test case 3). 
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Figure 39. Mach number contour predictions for the LeRC high area ratio nozzle using a viscous, multiple diffusing, 
chemically reacting analysis (Test case 4). 
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7. 1.4.2. Comparison of Experimental and Predicted Wall Pressures and Heat Fluxes 
In Figures 41 and 42 a comparison of the predicted (Test case 4) and measured static 
pressure and heat flux distribution in the supersonic region of the nozzle (the only region 
where measurements were made) is shown. The predictions and measurements are in 
excellent agreement. The ‘kink’ in the slope at the last three wall pressure measurements 
shown in Figure 41 is correctly predicted. 



Figure 41. Comparison of experimental wall static pressure measurements and 
predictions from the Tethys model. Test case 4. 

The heat flux profile shown in Figure 42 is observed to be in excellent agreement with 
experimental measurements, with a slight over prediction in the near throat region, fol- 
lowed by a slight under prediction, and then excellent agreement with experimental mea- 
surements as the exit of the nozzle is approached. 
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Figure 42. Comparison of experimental heat flux measurements and predictions from 
Tethys model. Test case 4. 

In Figure 43, an examination is made of the differences in wall heat transfer predic- 
tions between the three analysis techniques employed, using Test case 4 as the baseline 
heat flux case. Here it is observed that the difference between all of the analysis cases is 
rather small, being less than or about a maximum difference of 25 % between Test cases 2 
and 3 with Test case 4 tending to be within a 15 % difference of either of the other analysis 


cases. 
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Figure 43. Comparison of the difference in heat transfer predictions between the 
three different viscous flow analyses. 


7, 1,4.3. Prediction of Hydrogen Mm-Bnctiaa. Aipn g. ft e.Sfan 

In Figure 44 the predicted distribution of diatomic hydrogen mass fraction along the 
wall of the LeRC 1030:1 nozzle is presented for all of the viscous nozzle analysis cases 
performed. As expected, the concentration of hydrogen along the nozzle wall for the fro- 
zen analysis (Test case 2) is constant throughout the length of the nozzle. On the other 
hand, the mass fraction of hydrogen for the chemically reacting, ordinary diffusion case 
(Test case 3) increases slightly at the inlet due to the recombination of monatomic hydro- 
gen and then decays back towards the fully mixed, frozen mass fraction as the fluid under- 
goes a combination of species mixing and freezing of the chemical reaction mechanisms. 
The most pronounced difference in hydrogen mass fraction predictions is Test case 4, 
which considered Soret mass transfer. In this case it is seen that the hydrogen mass frac- 
tion immediately drops to 60 % of its inlet value, then slowly increases in mass fraction as 
the boundary layer thickens in the combustion chamber. As the nozzle begins to converge. 
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the hydrogen mass fraction continues to increase as the ceramic coating applied to the 
nozzle in this region serves to increase the wall temperature, with a corresponding 
decrease in the temperature gradient at the wall. The hydrogen mass fraction reaches a 
local maximum near the throat. Following this maximum, the hydrogen mass fraction rap- 
idly decreases as the Soret mass transfer term strengthens as a result of the rapidly 
decreasing wall temperature in this region. As the ceramic coated section ends, a discon- 
tinuous change in the rate of change of the mass fraction of the hydrogen is observed, fol- 
lowed shortly by a local minimum in hydrogen mass fraction. Although the wall 
temperatures continue to decrease beyond this region, the hydrogen mass fraction begins 
to increase as the temperature gradients at the wall begin to decrease and the boundary 
layer thickens, essentially remixing the fluid. This partial remixing process continues as 
the nozzle temperature gradient continues to decay throughout the extent of the expansion 
region of the nozzle. 



Figure 44. Hydrogen mass fraction predictions along the wall of the LeRC high 
area ratio nozzle. 
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7. 1 .4.4. Performance Predictions 

In Table 9, a comparison between the predicted performance parameters for the four 
test cases is shown. The procedure used to determine nozzle thrust and mass flow rate cor- 
respond as nearly as possible to that recommended in the JANNAF nozzle analysis proce- 
dure. A complete discussion of the JANNAF procedure can be found in Appendix M. The 
method involves the numerical integration of the pressure forces and the axial momentum 
flux across the throat of the nozzle, which in this case was taken to correspond to the first 
line that extends into the divergent portion of the nozzle, and to sum the pressure and drag 
forces along the divergent surface of the nozzle. 


Table 9. Tethys performance predictions. 



Test case 1 
(Frozen, 
Inviscid) 

Test case 2 
(Frozen, 
Viscous) 

Test case 3 
(Kinetics, 
Viscous) 

Test case 4 
(Kinetics, 
Viscous, 
with Soret 
and Dufour 
effects) 

Throat region thrust, lbf 

354.04 

351.97 

349.59 

349.51 

Throat region mass flow 
rate, m, lbm/sec 

1.1343 

1.1298 

1.1115 

1.1112 

Nozzle divergent region 
wall pressure thrust, lbf 

199.75 

203.86 

212.05 

212.20 

Nozzle divergent region 
shear forces, lbf 

0.00 

11.72 

12.96 

12.68 

Total vacuum thrust, Fy, lbf 

553.79 

544.11 

548.68 

549.03 

Specific impulse, F v /m, 
lbf-sec/lbm 

488.22 

481.60 

493.6 

494.1 

Thrust coefficient, 

f v /(p c a t ) 

1.959 

1.924 

1.941 

1.942 
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Focussing primarily on the final performance parameters of specific impulse and 
thrust coefficient, several conclusions can be drawn. The first is that finite rate kinetic 
effects seem to have the largest single effect on nozzle performance, with the performance 
difference between the frozen and kinetic analyses (Test cases 2 and 3, respectively) being 
about 2.5 %, with the frozen case being lower. Examination of Test cases 1 and 2 reveals 
that the viscous effects account for a 1.4 % decrease in nozzle performance. The effects of 
multidiffusion mass and heat transfer account for a predicted increase in nozzle perfor- 
mance of 0.1 %, as evidenced by comparison to Test cases 3 and 4. 

The mass flow rate predictions presented in Table 9 are also noteworthy. Examination 
of the predicted mass flow rates for Test cases 1 and 2 indicates that viscous effects result 
in a modest decrease in the mass flow rate through the nozzle. This is a boundary layer dis- 
placement effect. The primary difference in the predicted mass flow rates between Test 
cases 2 and 3 is that Test case 3 considered chemical kinetics in the analysis. Here it is 
seen that because of the release of sensible energy in the subsonic region of the nozzle in 
the chemically reacting analysis, the critical mass flow rate through the nozzle is about 
2 % less than in the chemically frozen analysis. As would be anticipated, the mass flow 
rate predictions for Test case 3 and 4 are essentially identical. 

7. 1 .4.5. Comparison of Experimental and Predicted Performance 

In Table 10, a comparison between the predictions of the Tethys fluid flow model. 
Test case 4, and experimental measurements is shown along with previous performance 
predictions that have been made with other models. The TDK prediction method (Refer- 
ence 1) involves an inviscid (one-dimensional/transonic/method of characteristics)- 
boundary layer type of analysis with chemical kinetics in the inviscid region while 
RPLUS (Reference 37) is a chemically reacting Navier-Stokes model which was executed 
using a 270X60 computational grid. 

A comparison of the predictions to the experimental measurements indicates that the 
Tethys flow model is in far better agreement with experimental measurements than the 
other computational models. The difference between the predicted and the experimentally 
measured thrust coefficient and specific impulse are approximately an order of magnitude 
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smaller in the Tethys model than the TDK or RPLUS nozzle prediction methods. 


Table 10. Experimental performance comparison of LeRC 1030:1 nozzle. 


Performance criteria 

Experimental 

measurements 

Tethys 
difference, % 

TDK 

difference, % 

RPLUS 
difference, % 

Thrust coefficient, 

f v /(p c a t ) 

1.946 

-0.21 

-2.5 

-1.4 

Specific impulse, 
F v /m, lbf-sec/lbm 

495.7 

-0.33 

-2.8 

-3.7 


Furthermore, since the experimental specific impulse reflects an inviscid, one-dimen- 
sional chemical equilibrium correction to account for incomplete combustion, the differ- 
ences between measured and predicted values of thrust coefficient and specific impulse 
should be nearly the same, with the only difference being the discharge coefficient of the 
nozzle, which will be near unity. As seen in Table 10, the proximity to which both the 
thrust coefficient and specific impulse are in agreement with experimental measurements 
lends additional credence to the accuracy of the Tethys model. 

















97 



A low thrust hydrogen/oxygen rocket engine was tested at the NASA LeRC as part of 
a nozzle contour evaluation program (References 38 and 39). The details of the nozzle 
geometry that will be used for comparative purposes is shown in Figures 45 to 48. In Fig- 
ure 45, a schematic of the nozzle geometry and test apparatus is shown, while in Figure 46 
a close-up view of the film cooling and injector apparatus is presented. In Figure 46 it is 
seen that a hydrogen film acts as a wall jet, the hydrogen being injected along the surface 
of the nozzle, and initially separated from the core by a solid walled sleeve. The sleeve 
from which the hydrogen film is injected consists of a channelled configuration with the 
channels merging into an axisymmetric configuration just prior to the end of the sleeve, as 
shown in Figure 47. A detail drawing of the throat region of the nozzle in shown in Figure 
48. 



Figure 45. 25 lbf film cooled rocket hardware and instrumentation schematic. 
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Spark Plug Fuel Splitting Washer 




Sleeve wall _ / 


0.015 in. 

0.05 in. 
0.03 in. 



30 Equally spaced channels 


3 - Sleeve inside diameter = 0.811 in. 

4 - Sleeve outside diameter = 0.870 in. 

5 - Channel outside diameter = 0.969 in. 

6 - Channel width = 0.055 in. 


Figure 47. Film cooled sleeve geometry. 
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Figure 48. Detail drawing of the 25 lbf film cooled nozzle. 


Similar to the LeRC 1030: 1 nozzle, this nozzle was tested in a near vacuum environ- 
ment. The rocket engine exhaust gasses flowed from the nozzle into a diffuser where the 
pressure was partially recovered. The final step in the removal of the exhaust gasses was 
accomplished by a combination of entrainment with ejectors and spray cooler condensa- 
tion (Figure 49). 



Figure 49. Experimental facility used for the 25 lbf film cooled nozzle tests. 
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The nominal operating conditions for the 25 lbf film cooled nozzle are given in Table 11. 


Table 11. Nominal operating conditions for the 25 lbf film cooled nozzle. 


Design parameter 

Nominal 

value 

- 

Chamber pressure, lbf/in. 

75.0 

Thrust, lbf 

25.0 

Overall mixture ratio 

8.0 

Core mixture ratio 

20.0 

Propellant mass flow rate, Ibm/sec 

0.08 


Experimental measurements of the following parameters were made: 

• Nozzle thrust 

• Fluid mass flow rate 

• Combustion chamber pressure 

• Nozzle wall static pressure 

• Nozzle wall temperature 

• Global heat transfer rate 

7.2.1. Analytical Assumptions for the 25 lbf Film Cooled Nozzle 
Three critical assumptions regarding the performance of the 25 lbf film cooled rocket 
were made so that analytical results could be generated in an expedient manner. The criti- 
cal assumptions made were that the flow was turbulent and could be reasonably modelled 
by an algebraic turbulence model, that the composition of the core inlet fluid flow of the 
nozzle was known, and that the radial and tangential velocity distribution at the nozzle 
inlet was zero for both the film and the core flow. 

The assumption of turbulent flow was based on the reasoning that the most crucial 
aspect of the analysis was the modelling of the film/core interaction in the combustion 
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chamber where the flow is almost certainly turbulent. On the other hand, the flow is also 
likely to experience a significant amount of laminarization, both along the wall and in the 
mixing layer, as it accelerates through the converging/diverging portion of the nozzle 
which would not likely be adequately modelled by an algebraic turbulence model. How- 
ever, since the nozzle is short, the overall shear stresses acting on the nozzle will be rela- 
tively minor and contribute only a small degradation of the total nozzle thrust. 

The assumption of a known chemical, thermal, and velocity distribution in the core 
inlet flow region was required. Because this was deemed to be a critical assumption in the 
analysis, two variations of the core flow inlet conditions were examined, one which 
assumed a fully mixed core and the other which was based on recent temperature and spe- 
cies measurements of a similar injector/combustor (Reference 40). 

Because the hydrogen film traversed through a three-dimensional channel configura- 
tion prior to injection/mixing with the core flow, it is certain that some distribution of 
radial and tangential velocity will exist at the inlet region of the hydrogen film cooled sec- 
tion. However, these components were assumed to be of relatively minor importance com- 
pared to the axial component of the velocity. Similarly, the oxygen and hydrogen 
propellants in the core region of the combustion chamber are injected radially and are also 
subjected to the effects of a blunt spark plug ignitor. Under these boundary conditions, the 
core flow will undoubtedly have some distribution of radial and tangential velocity for a 
distance downstream of the region where the propellants are injected. Unfortunately, ade- 
quate modelling of these conditions would require a three-dimensional analysis which is 
beyond the capabilities of the Tethys model. Instead, the sensitivity of the nozzle perfor- 
mance to the extent of core mixing was examined. 

7.2.2. Boundary Conditions for the 25 lbf Film Cooled Nozzle 

The boundary conditions applied to the flowfield simulation are presented in Table 
12. Several important and unique aspects of the boundary conditions employed are that the 
viscous, chemically reacting, reference plane characteristics method was employed at the 
inlet and exit boundaries with the exception of the solid wall region at the nozzle inlet. In 
this solid wall region, reference plane characteristics could not be used because the 
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turbulence model employed assumed that the distance from the wake producing body was 
relatively large. The characteristics technique employed along the nozzle centerline 
employed a chemically frozen, inviscid, waveline equation. It is also worthy of note that 
instead of the standard procedure of specifying the stagnation pressure at the nozzle inlet, 
the mass flux was specified, with the combustion chamber stagnation pressure that the 
nozzle was operating at being determined by the choking of the flow at the nozzle throat. 
Using this procedure of specifying the mass flow rates at the nozzle entrance, the flow split 
between the core and the film can be maintained exactly. 


Table 12. Boundary conditions applied to the 25 lbf film cooled nozzle. 


Inlet 

Exit 

Wall 

Centerline** 

Core 

Sleeve wall 

Film 

- 

- 

- 

MOC* 

^ = 0 
94 

MOC* 

MOC* 


MOC 

pu = 20.42 

u = 0 

pu = 4.0193 

MOC 

c 

II 

o 

* =0 
dTl 


v = 0 

v = 0 

MOC 

V = 0 

v = 0 

" - 9 <m 


h ° “ 582 lbm 

MOC 

T = T(x) 

*0 = 0 
dri 

Chemical 

equilibrium 

Ji = 0 

II 

© 

MOC 

Ji = o 

5(0 

5 — * = 0 

9ti 


* At subsonic contravarient 
velocities, ^exit ~ 0 * 15 psia. 

*At the points nearest the sleeve wall, the 

**Near the nozzle exit, reflective 

pressure gradient was assumed to be zero 

boundary conditions were 
applied. 
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The temperature distribution imposed on the nozzle wall was derived from interpola- 
tion of experimental data, where available, and is shown in Figure 50. The wall tempera- 
ture at the exhaust of the nozzle was not measured and was estimated by assuming that the 
temperature at the nozzle exit was equal to the water temperature that was used for cooling 
the nozzle walls. 



Figure 50. Temperature distribution used for the wall boundary condition of the 
Tethys simulation for the 25 Ibf film cooled nozzle. 

7.2.3. Computational Grid Used for the 25 Ibf Film Cooled Nozzle Study 
In Figure 51, the 91 X 51 computational grid used for the analysis of the 25 Ibf nozzle 
is presented. Since one of the fundamental aspects of the performance of this nozzle was 
the effectiveness of the hydrogen film coolant to mix and react with the oxygen rich core, 
a relatively large percentage of the computational grid points is contained within the sub- 
sonic portion of the nozzle. 



Radius, i 



0 . 1 . 2 . 3 . 4 . 5 . 

Axial distance,' in. 


Figure 51. Computational grid employed for the 25 lbf film cooled nozzle analysis. 
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The grid was generated by using the techniques described below. For 
1.2 < x < x Exit , 71 nozzle axial points were established using the following relation: 


(7.2.1) 


x = x 


Throat I 


i + 


sinh [x(£-B) ] 
sinh (xB) 


where B and x are as defined in Section 6.1, with the packing factor, x, used in the film 
cooled nozzle analysis being equal to 6.0. The total nozzle length and the distance to the 
throat are offset by 1.2 inches to accommodate packing in the combustion chamber section 
of the nozzle. For 0 < x < 1.2, the first 20 grid points were determined through the fol- 
lowing expression: 


(7.2.2) 


x 


1.2 


N- 1 
20 


1.5 


where N is the grid point number. The locations of the interior points were determined by 
extending a straight line from the centerline to the corresponding points established along 
the nozzle contour. The points along this line were packed near the nozzle wall using the 
following relation: 


(7.2.3) 


«> = 


r 

r nozzle 



1-T| 



+ 1 


where a=(3+l,b = (3-l, and (3 = 1.08, except near the throat region (from node num- 
ber 30 to node number 65), where the packing factor, (3, was linearly (in computational 
space) reduced to a value of 1.02 at the throat. 
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7.2.4. Analytical Study of the 25 lbf Film Cooled Nozzle 
Two different simulations were employed to examine the importance of stratified 
flow. One case assumed that the core of the flow was fully mixed and had an overall mix- 
ture ratio of 20. The other case assumed that the flow was highly stratified, as shown in 
Figure 52. The extent and profile of the core stratification was estimated by piecing 
together measurements on a similar injector (i.e., a square injector that operated at a core 
mixture ratio of 24, as discussed in Reference 40) which indicated that the centerline tem- 
perature of the core was about 3200 R and that the presence of hydrogen and oxygen was 
negligible near the sleeve wall. Applying these two pieces of information and further 
assuming that the temperature measurements were solely the result of incomplete mixing, 
a power law profile was generated. 



Figure 52. Oxygen to fuel mixture ratio profile assumed for the stratified core analysis. 

Two cases are presented in Table 13. Since the modelling of the film/core combustion 
process was deemed to be an important consideration, chemically reacting fluid flow, with 
eight chemical reactions, was common to all of the analyses. Turbulent flow was 
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considered for both of the cases presented herein. Since the flow was taken to be turbulent 
and the primary intent of this study was to examine the mixing of the fluid core and the 
film, all mass diffusion terms and heat transfer modes were employed. Likewise, the tur- 
bulent Prandtl number used in this analysis corresponds to experimental measurements 
obtained in a free turbulent shear flow environment (Reference 41). 


Table 13. Analysis methods used for the 25 lbf film cooled nozzle. 


Tethys 

Test case 1 

Test case 2 

analysis 

(fully mixed core) 

(stratified core) 

Number of species considered 

6 (H 2 0, H 2 , 0 2 ,OH, 

6 (H 2 0, H 2 , 0 2 , OH, 


H, 0) 

H, O) 

Thermal diffusion 

Yes 

Yes 

Interdiffusion heat transfer 

Yes 

Yes 

Dufour heat transfer 

Yes 

Yes 

Kinetics 

Yes 

Yes 

Number of chemical reactions 

8 

8 

Turbulence 

Yes 

Yes 

Artificial diffusion 

Nonlinear, second 

Nonlinear, second 


order model, 100 % of 

order model, 100 % of 


recommended value. 

recommended value. 

Artificial diffusion added to 

Inlet- Yes 

Inlet-No 

MOC equations 

Exit- Yes 

Exit- Yes 


Centerline- Yes 

Centerline- Yes 

Turbulent Prandtl number 

0.5405 

0.5405 

Turbulent Lewis number 

1.0 

1.0 


The algebraic turbulence model used in the film cooled nozzle simulation employed 
the Baldwin-Lomax inner region model and a free jet turbulence model in the outer 
region. The free jet turbulence model used in this analysis was one applied by Gortler 
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(Reference 42) for the derivation and verification of the velocity profile of a mixing layer: 
(7.2.1) (i t = 0.0014p (x + Ax) Au 

where the velocity defect, Au , applied in this analysis corresponded, in an approximate 
manner, to the average velocity difference between the core and the film coolant at the 
inlet For consistency, the velocity defect used for both cases was the same, 500 ft/sec. The 
distance, x, in the above correlation corresponded to the physical distance of a given point 
to the inlet, along a constant T) line. The critical length calculation is offset by an addi- 
tional distance. Ax, of 0.144 inches (0.012 feet). The value of this offset distance corre- 
sponds to the theoretical distance at which the thickness of the mixing layer would be 
equal to the thickness of the sleeve wall, as illustrated in Figure 53. 


Film coolant flow 





Figure 53. Mixing layer origination point. 

While very simple in nature, the above free-stream turbulent mixing layer correlation 
captures the fundamental features of the mixing process. The ‘freezing’ of the turbulence, 
due to flow acceleration in the convergent-divergent portion of the nozzle (see Reference 
43) is also reasonably well modelled since the velocity defect is a constant, while the den- 
sity decreases appreciably (about two orders of magnitude) as the flow accelerates through 
the nozzle. It is worth noting that the standard free-stream turbulence model employed in 
the Proteus model yields an identical expression for the turbulent viscosity calculation of a 
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mixing layer (if the theoretical maximum shear stress value is inserted), the only differ- 
ence being that the constant is 40 % lower than applied in this analysis and the velocity 
defect is based on conditions along a constant £, line. 

7.2.5. Analytical Results and Comparison with Experimental Measurements 
To examine and compare the analytical predictions with experimental measurements, 
an examination of the Mach number contours, fluid velocity vectors, species concentra- 
tions, pressure distribution, and wall heat flux predictions were made on a local scale. The 
global heat transfer rate and nozzle performance were compared for the two test cases and 
measured data. 

7.2.5. 1. Mach Number Contour Predictions 

Due to the conical nature of the nozzle divergence, one would expect that modest 
shock waves would exist in the supersonic region of the nozzle. In Figures 54 and 55, the 
predicted Mach number contours for the fully mixed and stratified core flow test cases, 
respectively, verify what would be intuitively assumed to exist. A rather weak shock wave 
coalesces at the nozzle centerline and causes the Mach number to decrease from a value of 
around 3.0 to a value somewhat greater than a Mach number of 2.5. The Mach number 
contours are essentially insensitive to the nature of the flow stratification, appearing 
almost identical in both the fully mixed and the stratified core cases. The Mach number 
contours are highly distorted in the high supersonic region of the nozzle, being much 
higher (at a given axial location) near the wall than near the centerline. It is evident that 
the shock waves had a far greater effect near the nozzle centerline region than near the 
wall region. 




Figure 54. Mach number contour predictions for the 25 lbf film cooled nozzle using the fully mixed core assumption. 


o 



Figure 55. Mach number contour predictions for the 25 lbf film cooled nozzle using the stratified core assumption. 
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1.2.52. Fluid Velocity Vectors 

The complicated nature of the film cooled nozzle system and the conical geometry of 
the supersonic portion of the nozzle might be expected to yield rather complex velocity 
vectors. However, the predictions from the stratified flow analysis, illustrated in Figures 
56 and 57 with every fifth r\ line and every third £, line printed, indicate that the flow is 
very regular in nature, even at the boundaries. 




Figure 56. Velocity vector predictions for the 25 lbf film cooled nozzle using the stratified core assumption. 



Figure 57. Throat region velocity vector predictions for the 25 lbf film cooled nozzle using the stratified core assumption. 
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7.2. 5.3. Predicted Steam Mass Fraction Contours 

In Figures 58 to 61, a contour mapping of the mass fraction of steam (H 2 O) is shown 
for the fully mixed core analysis and the stratified core analysis. In examining these fig- 
ures, numerous flow features become very evident The first is that in both the fully mixed 
and the stratified core cases, fluid mixing is essentially nonexistent in the supersonic 
region of the nozzle. This is entirely as expected, since the fluid speeds increase rapidly in 
the transonic portion of the nozzle, while the turbulence quantities tend to ‘laminarize’. 
Another important observation is that the stratified flow case exhibits a much lower mass 
fraction of steam in the core region of the flow throughout the nozzle. Since the creation of 
steam is associated with a large increase in thermal energy, one can expect a significantly 
higher performance from the fully mixed core case. As shown in Section 7 .2.5.6 , this is 
exactly the situation. 





0.4 


Figure 59. Steam mass fraction contours in the throat region of the 25 lbf film cooled nozzle using the 
fully mixed core assumption. 
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7.2.5.4. Predicted and Measured Pressure Profiles 

Analytical predictions of nozzle wall pressure and experimental pressure measure- 
ments are presented in Figures 62 and 63 for the fully mixed core analysis and the strati- 
fied core analysis cases, respectively. In both of these cases, the wall pressure predictions 
correlate extremely well with measured data. However, in both cases the predicted pres- 
sures are slightly higher than the experimental measurements, with the fully mixed core 
case being somewhat higher than that predicted in the stratified core case. Most of the dis- 
crepancy between predicted and measured wall pressures can be attributed to the fact that 
the predicted combustion performance was somewhat better (from 6 % to 11 %, as dis- 
cussed in Section 7.2.5.6) than occurred in the experimental program. As the mass flow 
rates were used as an input to this analysis, the higher combustion performance is reflected 
in a higher stagnation pressure at the nozzle inlet. 

In both cases, a significant discontinuity in the slope of the pressure distribution is 
predicted at the location where the diverging nozzle cone attaches to the curved transonic 
portion of the nozzle. The crispness in the modelling of this effect is noteworthy and 
would be fully expected to exist in actual test conditions, although the coarseness of the 
pressure measurements gave no indication of the presence of this effect. 



Figure 62. Comparison of experimental and predicted wall pressure distributions in 
the 25 lbf film cooled nozzle using the fully mixed core assumption. 
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Figure 63. Comparison of experimental and predicted wall pressure distributions in 
the 25 lbf film cooled nozzle using the stratified core assumption. 


12 . 5 . 5 . Predicted and Experimental Heat Transfer Rates 

Global measurements of nozzle heat transfer and local measurements of nozzle wall 
temperature were made. The predicted nozzle wall heat fluxes are shown in Figures 64 and 
65 for the fully mixed core and stratified core cases, respectively. While a high level of 
confidence in the quantitative accuracy of the predictions cannot be expected due to a vari- 
ety of reasons (which are discussed on page 124), numerous aspects of the heat flux profile 
predictions are of importance. 

An extremely interesting observation is that the peak in the nozzle heat fluxes is 
appreciably upstream of the nozzle throat, with the peak heat flux location being about 0.2 
inches upstream of the nozzle throat. This location is well upstream of the small amount 
that is normally experienced due to nozzle curvature which causes the sonic line to curve 
upstream of the throat. This maximum heat flux point would be expected to occur at about 
0.04 inches upstream of the nozzle throat (Reference 23, page 100). Instead, the peaking 
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of the heat flux profile represents the prediction of injector streaking. The streaking occurs 
because the momentum of the hot film core is significantly larger than that of the hydrogen 
film coolant, being about 20 times greater, and the core flow has the tendency towards 
direct impingement onto the wall since the low momentum of the coolant has little ability 
to force the core flow to turn with the geometry. Instead, the fluid streamlines become con- 
centrated very near the wall of the nozzle. 



Figure 64. Wall heat flux predictions for the 25 lbf film cooled nozzle using the fully 
mixed core assumption. 


123 



Figure 65. Wall heat flux predictions for the 25 Ibf film cooled nozzle using the 
stratified core assumption. 

Another interesting aspect of the heat flux predictions shown in Figures 64 and 65 is 
the local minimum, followed shortly by a local maximum that occurs just downstream of 
the throat In fact the sharp increase in heat flux corresponds precisely to the region where 
the geometric transition from the circular arc throat to the conical divergence occurs. 

The streaking aspect of the injector is used in Figure 66 to present the wall tempera- 
ture distribution which would likely exist in the region just upstream of the throat. These 
predictions (for the stratified core case and the fully mixed core case) are compared to the 
straight-line-fit temperature distribution as published in a report (Reference 38) detailing 
the experimental results. Here it is seen that the actual wall temperatures may have been in 
excess of 100 degrees Rankine hotter than suggested with a straight line fit of the data. 
Rather surprisingly, the temperature increase is mostly independent of the predictive tech- 
nique employed. This is the result of scaling the coolant channel resistances to fit the 
experimentally measured wall temperatures. Hence, the temperature increase is a 
reflection of the relative increase of the heat flux in this region and is independent of its 
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actual magnitude. 



Distance from throat, in. 

Figure 66. Streaking effects on the wall temperature distribution of the 25 lbf film 
cooled nozzle. 

The global heat transfer rate measurements are presented in Table 14. Surprisingly, 
the stratified core analysis predicts heat transfer rates that are within 2.6 % of experimental 
measurements. Because the turbulence model was rather simplistic in nature and the tur- 
bulent Prandtl number was tailored to accurately reflect free-stream mixing conditions, in 
addition to a rather coarse, loosely packed grid and the fact that the constituency of the 
stratified core itself is only modestly well defined, it is most likely that this high level of 
agreement between analysis and experiment is a fortuitous result. The predicted heat 
transfer rates of the fully mixed core nozzle are substantially higher than experimental 
measurements or that predicted in the stratified core analysis. Undoubtedly, this is the 
result of the much higher heat flux peak at the region upstream of the throat which in turn 
is likely due to the existence of unreacted oxygen in close proximity to hydrogen. 



125 


Table 14. Comparison of measured and predicted global heat transfer rates. 



Total heat transfer, 
Btu/sec 

Experimental measurements 

23.6 

Stratified core analysis 

23.0 

Fully mixed core analysis 

31.5 


12 . 5 . 6 . Predicted and Measured Performance Parameters 

One of the most important requirements of a nozzle analysis is the ability of the 
model to adequately predict nozzle performance parameters such as specific impulse and 
thrust coefficient. In the case of this particular film cooled nozzle analysis, this was espe- 
cially difficult since the chemical composition of the core flow was only modestly well 
defined. Further, the mixing effects between the film and the fluid core were essential to 
predict specific impulse and characteristic velocity and would seem to depend heavily on 
the nature of the turbulence model. In light of these seemingly difficult conditions, the crit- 
ical performance parameters regarding nozzle performance are compared in Figures 67 
and 68 and in Table 15 between the two nozzle analysis cases and experimentally mea- 
sured results, at a mixture ratio of 8.0. 

The measured specific impulse for two slightly different injectors (the difference 
being a slight difference in the size of the gap between the spark plug base and the com- 
bustor walls) used in the test program and the predictions from the Tethys code are shown 
in Figure 67. Here it is seen that both of the Tethys analysis cases tend to overpredict the 
performance of the nozzle, but that the performance predictions of the stratified flow case 
are substantially closer to experimental results than the fully mixed core case. Since the 
sensitivity of the nozzle to slight changes in the injector geometry caused such a large dif- 
ference in the experimentally measured performance, being about a 10 lbf-sec/lbm differ- 
ence at a mixture ratio of 8.00, it is highly probable that the stratified core flowfield 
assumed at the inlet of the Tethys analysis underestimated the extent of flow stratification. 




Vacuum thrust coefficient 
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Figure 67. Comparison of experimental and predicted specific impulse for the 25 lbf 
film cooled nozzle study. 



Figure 68. Comparison of experimental and predicted thrust coefficient for the 25 lbf 
film cooled nozzle study. 
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A comparison of the measured and predicted thrust coefficients for the film cooled 
nozzle is shown in Figure 68. Here it is seen that the Tethys predictions for either the strat- 
ified core or the fully mixed core are well within the experimental disparity observed 
between the two injectors that were tested. 

A final observation can be drawn by relating the trends of Figure 67 and those of Fig- 
ure 68, always focussing on the results obtained near the mixture ratio of 8.0. In Figure 67, 
it is observed that the experimentally measured specific impulse is about 3 % higher for 
injector SN. 02 as compared to the specific impulse measured for injector SN. 03. In Fig- 
ure 68, the thrust coefficient observed for injector SN. 02 appears to be about 1 % higher 
than the thrust coefficient of injector SN. 03. These observations are in contrast to those of 
the Tethys code results where it is observed that a 5 % difference in specific impulse 
between the fully mixed core assumption and the stratified core assumption relates to a 
predicted difference in vacuum thrust coefficient of only 0.2 %. While the experimental 
and predicted trends of higher specific impulse corresponding to higher thrust coefficients 
is in agreement, the disparity between the extent to which this occurs suggests that the 
flow is not only more stratified than used in the Tethys analysis, but that the qualitative 
nature of the profile is also inaccurate. 

In Table 15, a precise comparison between the measurements obtained using the 
lower performing injector (i.e., injector SN. 02) and the prediction of the Tethys model is 
made. Here it is seen that the specific impulse and characteristic velocity predictions of the 
stratified flow case are about 6 % higher than was experimentally measured while the 
thrust coefficient prediction is within 0.5 % of the experimentally measured value. 

In general, the difference between experimental measurements and the stratified core 
predictions regarding engine performance are about equal to one-half the difference 
between experimental measurements and the fully mixed core predictions. The exception 
to this statement is the difference in experimental measurements concerning the nozzle 
thrust coefficient, where the difference between measured and predicted values in the 
stratified core analysis is about 50 % closer to experimental measurements than the fully 
mixed core analysis. 
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Table 15. Comparison of measured and predicted performance parameters. 


Engine 

performance 

parameter 

Experiment 

Fully mixed 
core 
analysis 

Stratified 

core 

analysis 

Difference 
(Analysis / 
Experiment), % 

Specific impulse, 
F v /m, lbf-sec/Ibm 

304.1 

339.9 

323.5 

Mixed: + 11.5 
Stratified: + 6.4 

Characteristic veloc- 
ity, P c A T /m, ft/sec 

5947. 

6599. 

6296. 

Mixed: + 11.0 
Stratified: + 5.9 

Vacuum thrust coeffi- 
cient, F v / (P C A X ) 

1.645 

1.657 

1.653 

Mixed: + 0.73 
Stratified: + 0.49 
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7.3. Transpiration Cooled Rocket Nozzle 

A hydrogen/oxygen transpiration cooled rocket engine was tested at the NASA LeRC 
in the 1970’s (Reference 44). The purpose of the program was to verify the ability of a 
hydrogen coolant to significantly reduce heat fluxes in the throat region of a rocket nozzle 
and to demonstrate that a new technique of creating a transpired region by a series of fine 
platelet stacks produced a reliable, efficient, and economically viable means of transpira- 
tion cooling rocket throat sections. Figures of the baseline rocket engine and the transpira- 
tion cooled spoolpiece section are shown in Figures 69 and 70. As seen in Figure 69, the 
demonstration tests for hydrogen transpiration cooling were not performed on a conven- 
tional converging-diverging rocket nozzle. Instead tests were performed with a plug and 
spool engine, with the transpiration coolant being injected into the main stream of the 
rocket exhaust gasses only along the spoolpiece and limited to a region 0.5 inches in 
length (see Figure 70), the region ending at the throat. The overall length of the spoolpiece 
section was 4 inches, and the gap between the spoolpiece and the plug was 0.25 inches at 
the throat and 0.5 inches in the combustion chamber. The inside diameter of the spool- 
piece was 2.6 inches. 



E 



Figure 69. Plug and spool rocket engine tested at NASA LeRC. 
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Figure 70. Schematic of the transpiration cooled spoolpiece section. 

The nominal operating conditions of the nozzle are given in Table 16. The mixture 
ratio and chamber pressure were nearly constant throughout the test series, while the 
hydrogen transpiration mass flow rates were varied from 0.045 lbm/sec to 0.075 lbm/sec. 


Table 16. Nominal operating conditions of the transpiration cooled rocket nozzle. 


Design Parameter 

Nominal 

Value 

Chamber pressure, lbf/in. 

600.0 

Oxygen / Fuel mixture ratio 

6.0 

Transpired hydrogen mass flow rate, lbm/sec 

0.045 

Transpired section wall temperature, R 

900.0 
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7.3.1. Analytical Modeling of the Transpiration Cooled Rocket Nozzle 
Several important assumptions were made concerning the ultimate predictive capa- 
bilities of the Tethys model in predicting the flowfield of the transpiration cooled rocket 
nozzle. Due to the relatively incomplete nature of the experimental data reporting tech- 
niques and the significant computational advantages, the hydrogen/oxygen propellant sys- 
tem was assumed to be a frozen mixture of steam and excess hydrogen. Furthermore, the 
temperature of the transpired coolant was presumed to be equal to the wall temperature of 
the copper platelet stack through which the fluid flowed. While no direct measurements of 
the fluid temperature being injected were made, a heat transfer analysis performed during 
the experimental program indicated that the wall temperatures might be as much as 300 
degrees Rankine hotter than the coolant fluid. Instead, the estimated temperature of the 
coolant fluid was used to represent the wall temperature in the transpired section, reason- 
ing that the enthalpy difference of the coolant is far more critical than the slightly different 
conduction heat transfer rates that may occur as the result of the wall temperature being 
suppressed. 

7.3.2. Boundary Conditions Employed for the Transpiration Cooled Rocket Nozzle 
The boundary conditions applied in the Tethys model simulation are presented in 
Table 17. As shown in the table, the reference plane characteristic method (MOC) was 
used at all of the boundaries. In this transpiration cooled analysis, the wall temperature 
distribution was specified throughout the transpired and solid wall regions. Hence, the 
mass flow rate requirements needed to yield the prescribed wall temperature was a vari- 
able to be determined. 
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Table 17. Boundary conditions applied to the transpiration cooled rocket nozzle. 


Inlet 

Exit 

Plug Wall 

Spool Wall 

Solid 

Transpired 

P 0 = 600 lbf/in. 2 

MOC* 

MOC 

MOC 

MOC 

MOC 

MOC 

* =0 
3r| 

u = 0 

u = 0 

v = 0 

MOC 

o 

II 

c 

> 

V = 0 

pvAh + q" = 0 

ST 

o 

II 

1 

t— » 

o^W 

3 e 

MOC 

^° = o 

T = T(x) 

H 

II 

H 

<o H2 = 0.035714 

MOC 

3c0tt 

^ 2 = 0 

o 

II 

X 

" — > 

Jh 2 + P VC0 h 2 = P v 

♦At subsonic contravarient velocities, the exit 
pressure was set equal to 61.5 lbf/in. 2 




An additional consideration that required investigation was the effect of the tran- 
spired surface wall on the flow stream. Due to excessive time in an electropolishing bath, 
the plates employed for the construction of the transpiration region became rounded, 
resulting in an extremely rough surface, as illustrated in Figure 71. Consequently, the 
effect of this roughened surface was examined by conducting a smooth walled coolant 
section analysis and a rough walled coolant section analysis. 

The extent of the rounding caused by the excessive polishing process was measured 
to be 12 RMS on “Taly-Surf ’ measurement equipment and corresponded very well with 
magnified photographs. 
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Rocket engine exhaust gasses 



Transpiration cooled plates, as designed 


Rocket engine exhaust gasses 



Transpiration cooled plates, as tested 


Figure 71. Schematic of surface roughening caused by excessive electropolishing of the 
transpiration cooled rocket nozzle. 


Table 18 summarizes the different assumptions and methodologies employed for the 
two analyses. 


Table 18. Analysis methods used for the transpiration cooled rocket nozzle. 


Analysis Parameter 

Smooth Wall 

Rough Wall 

Soret mass transfer 

Yes 

Yes 

Dufour heat transfer 

Yes 

Yes 

Turbulence model 

Cebeci-Smith 

Cebeci-Smith 

Turbulent Prandtl number 

0.9 

0.9 

Turbulent Lewis number 

1.0 

1.0 

Surface roughness, transpired region 

None 

12 mils (0.001 ft) 

Artificial diffusion (main flowfield and 

Second-order non- 

Second-order non- 

MOC equations, except MOC equations 

linear model: 2X 

linear model: 2X 

applied along the spoolpiece wall) 

recommended level 

recommended level 
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The spoolpiece temperature distribution used for the transpiration cooled nozzle 
study is shown in Figure 72. In the solid walled region, the temperature was established 
from experimental measurements. In the transpiration cooled region, a combination of 
experimental measurements and temperature correction factors were used to estimate the 
coolant temperature distribution. The reason for the depression of temperature in the cen- 
ter of the cooled region is possibly due to the nature of the hardware/engine design system 
which, when heated, preferentially closes those coolant slots that are located nearest the 
solid wall region. Consequently, the coolant flow will preferentially escape through the 
center of the transpiration coolant region. 



Figure 72. Wall temperature used for analysis of the transpiration cooled 
rocket nozzle. 


7.3.3. Computational Grid Used for Analysis of the Transpiration Cooled Rocket Nozzle 
A computational grid with 51 axial and 31 radial grid points was used for the analysis 
of the transpiration cooled rocket nozzle. The distribution of grid points is presented in 
Figure 73, where it is observed that a significant amount of axial packing was employed in 
the region enclosing the transpiration cooled section, and the grid was heavily packed near 
the no-slip spoolpiece surface. 
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Figure 73. Computational grid employed for the analysis of the transpiration cooled rocket nozzle. 
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The computational grid was established using algebraic relations. To establish the 
axial location of the computational lines, the grid was constructed by first assuming that 
only 31 axial lines were to be constructed, with the following relation: 


(7.3.1) 


x = x 


Throat 


l + 


sinh [ x(^-B) ] 
sinh (xB) 


where B and X are as defined in Section 7.1, with the packing factor, x, equal to 1.001, 
which, in essence, resulted in an evenly distributed grid. In the region encompassing the 
transpiration cooled section, the axial grid point spacing was 20 % of the spacing used in 
the inlet and exit regions. The extremely tight radial packing employed in this analysis 
was accomplished using the following relation: 


(7.3.2) 


<t> = 


a - 



1 — T| 


1 + 



1-T1 


with <j> being equal to the ratio of the distance from the plug wall divided by the distance 
from the plug wall to the spoolpiece. The constants a and b given in Equation (7.3.2) are 
equal to P +1 and P — 1, respectively, with the packing factor, p, being equal to 1.0004 in 
the smooth surface regions and equal to 1.0008 in the rough surface regions. The relax- 
ation of the packing factor, p, for the rough walled region in the rough walled analysis 
reflects the fact that significantly less packing is required for a rough walled analysis, since 
the viscous sublayer does not exist. 

Near the wall of the plug surface, enhancement of the radial packing was performed 
by overwriting the grid points generated in Equation (7.3.2) with the correlation that is 
presented below. The enhancement was applied to the eight points nearest the surface of 
the plug wall, which corresponds to the T] = 1 surface. 




r-r 22 


N - 22 


1.0 -0.5 


N - 22 
31-22 


(7.3.3) 
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7.3.4. Analytic Results and Comparison with Experimental Measurements 

Since the primary intent of this study was to examine the capabilities of the Tethys 
model for establishing the mass flow rates required to maintain a rocket engine at a pre- 
scribed temperature, the results to be presented will focus primarily on predictions and 
measurements made in the transpiration cooled region, although some levels of global pre- 
dictions will be shown to establish the consistency of the Tethys model. Specifically, in the 
transpired region, the balance in heat flux and mass flow rate will be shown for the two 
different analyses and will be compared to the overall experimental measurements. Addi- 
tional results pertaining to predictions involving the flowfield Mach number contours and 
surface parameters such as heat fluxes and hydrogen mass fractions will also be shown. 

7.3.4. 1 . Mach Number Contour Predictions 

The extremely gentle curvature of the subsonic and supersonic regions involved in 
this geometry would be expected to yield predictions that are nearly one-dimensional in 
nature. The predicted Mach number contours, as illustrated in Figures 74 and 75, support 
this reasoning. The Mach number contours have only a modest amount of radial curvature 
associated with them. As would also be anticipated, the smooth and rough walled analyses 
yield essentially identical Mach number contours, although extremely close examination 
of the contours indicates that the subsonic Mach number contours of the rough walled 
analysis are more tightly packed than those of the smooth walled analysis. This is the 
result of the increased mass addition experienced with the rough walled analysis. In turn, 
the increased mass flow rate of the rough wall displaces the boundary layer more than the 
smooth walled counterpart and consequently, the effective area ratio available for flow 
acceleration is sightly altered. 






Figure 75. Predicted Mach number contours for the transpiration cooled plug and spool rocket nozzle using a 
rough wall analysis. 
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7.3. 4.2. Hydrogen Wall Mass Fraction Prediction 

The mass fraction of hydrogen along the wall of the spoolpiece is a parameter of 
interest because of the relation of this term to mass fluxes and the significance of the Soret 
mass transfer term as compared to other diffusion terms. It is not entirely clear a priori 
what the effect of turbulence is on the Soret mass diffusion term. On the one hand, the tur- 
bulence tends to drive the species towards a fully mixed condition, while on the other hand 
turbulence creates extremely large temperature gradients very near the wall, thereby mak- 
ing the Soret term also large. In Figures 76 and 77, the predictions for the mass fraction 
distribution of hydrogen are shown for the smooth and rough transpiration surface cases, 
respectively. The profiles look very similar in nature, with the exception being that the 
hydrogen mass fraction in the rough walled analysis is significantly higher than that of the 
smooth walled case in the region downstream of the transpiration section. Undoubtedly, 
this is the result of the additional hydrogen required as a coolant in the rough wall analy- 
sis. 

A final observation is that upstream of the transpiration cooled section, the mass frac- 
tion of hydrogen on the wall is actually about one-half that of the core, which is a clear 
indication of the strength of the Soret mass transfer term even in the presence of turbulent 
flow. 



■2 
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Figure 76. Hydrogen wall mass fraction predictions using a smooth walled analysis. 



Figure 77. Hydrogen wall mass fraction predictions using a rough walled analysis. 
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7.3.4.3. Transpiration Coolant Mass Flux Predictions 

The critical criteria upon which the success of the transpiration cooled nozzle analy- 
sis will be based is the ability to adequately predict the mass flux required to maintain a 
wall at a prescribed temperature. Because the mass flow that enters a transpired region 
may exit by diffusion (either ordinary or thermal) and advection, the relative importance 
of each term is of interest. 

In Figure 78, the predictions of mass flux in the transpiration cooled section of the 
nozzle are presented for the two cases performed along with the averaged mass flux mea- 
sured in the experimental program. Numerous observations become readily apparent One 
of the most obvious is that the mass flux required of a smooth transpiration cooled wall are 
appreciably lower than was experimentally observed. On the other hand, the rough walled 
analysis appears to track the experimental mass flux rates rather well until reaching the 
center of the transpiration section, where the prescribed wall temperature begins to 
increase. At this point, the mass injected into the fluid stream acts as a film coolant and 
serves to reduce transpiration requirements until reaching the edge of the transpiration 
cooled region where the highly accelerating flow of the rocket engine serves to dramati- 
cally increase coolant requirements, although the extent to which this occurs may be over- 
exaggerated because of the discontinuous nature of the adjoining boundary conditions. 
Similar trends are observed in the smooth walled analysis but to a much less extreme. 
Table 19 summarizes the differences between the observed and predicted transpiration 
mass flow rates in the plug and spool rocket engine. Here it is observed that the smooth 
walled prediction is nearly 50 % that of the experimental measurements, while the rough 
walled prediction is much closer to experimental results, being only 13 % lower than the 
experimentally measured mass flow rate. 
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Table 19. Total coolant flow rate requirements. 



Transpiration wall analysis 

Experiment 

Difference 


Smooth 

Rough 

— 

(Analysis/ 
Experiment), % 

Hydrogen coolant mass 
flow rate, lbm/sec 

0.022 

0.039 

0.045 

Smooth: -51.0 
Rough: -13.0 


The individual contributions of mass flux into the fluid stream are shown in Figure 79 
and Figure 80 for the smooth and rough walled analysis cases, respectively. 
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Figure 79. Predicted hydrogen mass fluxes along the wall of the transpiration 
cooled plug and spool rocket engine using a smooth wall analysis. 

Examination of the smooth walled analysis case yields many interesting features. The 
first is the importance of the Soret mass flux term in both the solid walled and transpiration 
cooled sections of the nozzle. As the term is always positive, the Soret term drives the 
lighter fluid (hydrogen) away from the wall and towards the mainstream of the rocket 
engine exhaust gasses. In the solid walled regions of the flow, the positive Soret hydrogen 
flux term is balanced by the ordinary diffusion terms to create a no net flux region. On the 
other hand, in the transpiration cooled section of the engine, the ordinary diffusion flux is 
directed in the positive direction with the net result being that all hydrogen mass flux 
terms, advective and diffusive, act to drive the hydrogen into the flow stream. 

The rough wall results exhibit most of the same trends as the smooth wall case with 
the exception that in the transpiration cooled section of the engine the contributions of 
hydrogen mass flux as the result of ordinary and thermal diffusion are very small com- 
pared to the contributions of the advective and turbulent diffusion terms, which would be 
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expected because of the extremely high level of surface roughness used in the analysis. 



Figure 80. Predicted hydrogen mass fluxes along the wall of the transpiration 
cooled plug and spool rocket nozzle using a rough wall analysis. 


7.3.4. 4. Heat flux predictions 

The individual contributions of the diffusion heat transfer components are presented 
in Figures 81 and 82 for the smooth and rough walled cases, respectively. The heat trans- 
fer rates due to conduction are seen to vary dramatically throughout the nozzle length. In 
the smooth wall results, the conduction flux decays significantly in the transpiration 
cooled section of the nozzle, being about 20 Btu/(in. -sec) just prior to the transpiration 
cooled section and decaying to values of less than 10 Btu/(in. -sec) in the transpiration 
cooled section. The effect of this film cooling is actually even greater than suggested as 
the conduction heat flux would normally (i.e., if the wall were solid) be expected to 
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increase by 50 % in this region as the result of flow acceleration (Reference 45). In the 
rough walled analysis, the conduction heat flux is seen to oscillate much more signifi- 
cantly with the general tendency of the conduction flux to rise rapidly at the initial section 
of the roughened/transpired section, proceeded by a sharp decay in temperature gradients 
in the final half of the transpiration cooled section as the prescribed wall temperature 
begins to rise. 



Figure 81. Predicted heat flux into the mainstream of the plug and spool rocket 
nozzle using a smooth wall analysis. 
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Figure 82. Predicted heat flux into the mainstream of the plug and spool rocket 
nozzle using a rough wall analysis. 

Unlike conduction, the flux of heat as the result of interdiffusional heat transfer is 
observed to be directed into the main stream of the flow. Since the expression for interdif- 
fusional heat transfer is the product of the mass flux of hydrogen (which is positive in the 
transpiration cooled section) and the difference between the enthalpy of hydrogen and 
steam, the diffusion heat flux will be rather large as the enthalpy of hydrogen is near zero 
at the prescribed spoolpiece wall temperatures but the enthalpy of steam is very large and 
very negative at the same wall conditions. 

A final observation to be made is that the flux of heat as the result of Dufour energy 
transfer is negligible in both of the analyses performed. 
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SECTION 8 
CONCLUSIONS 

A model for the prediction of rocket nozzle performance and heat transfer has been 
developed and demonstrated for a rather wide spectrum of nozzle geometries and bound- 
ary conditions, with the overall conclusion to be drawn that the flowfield model is both 
extremely versatile and extremely accurate for propulsive nozzle flow simulations. In 
many aspects, the model is the most comprehensive in the world, as the flow modelling is 
conservative for both inviscid and viscous flows, and contains a comprehensive capability 
to model multispecies ordinary diffusion, thermal diffusion, Dufour heat transfer, and 
chemical kinetics, and is further enhanced by the use of a reference plane method of char- 
acteristics at computational boundaries. 

The Tethys methodology consistently predicted rocket nozzle thrust coefficient to 
values within 0.5 % of experimental measurements, and it predicted heat fluxes extremely 
well under a variety of different nozzle geometries and operating conditions. Furthermore, 
wall static pressure predictions agreed extremely well with experimental measurements, 
with the Tethys model predicting kinks and jumps in wall pressure that were either experi- 
mentally verified or fully anticipated due to theoretical considerations. 

Additional conclusions to be drawn from this modelling effort as regards the model- 
ling of chemical rocket engines in the demonstration studies are contained below. 

Soret mass transfer considerations can significantly affect species mass transfer and 
mass concentrations. In a solid wall region, the thermal mass transfer term resulted in sig- 
nificant changes in species concentrations predictions while yielding only a slight change 
in heat flux predictions. This is reasoned to be due to the fact that the species flux is zero 
near a solid wall region. Consequently the potential impact of Soret diffusion on heat 
transfer is by the change in fluid transport properties, and in the one demonstration case 



150 


examined, the change was not appreciable. In a transpiration cooled rocket nozzle, the 
Soret mass transfer term was a significant diffusion parameter at the wall, provided the 
surface of the rocket was smooth. 

Dufour heat transfer is relatively small. Since the diffusion coefficient in interdiffu- 
sional heat transfer contains the difference in enthalpy between two species while that of 
Dufour heat transfer is the difference in molecular weight, a system comprised primarily 
of excess fuel and combustion products will likely have a far greater interdiffusional heat 
transfer rate than Dufour heat transfer. 

Chemical kinetics has a significant effect on nozzle performance. In the 1030:1 high 
area ratio nozzle examined, the performance difference between finite rate chemical kinet- 
ics and frozen flow was significantly greater than for other effects. 

The viscous, chemically reacting, reference plane method of characteristics is very 
well suited for chemically reacting rocket nozzle inlets and exhausts. Even in the viscous 
regions, the reference plane characteristics method performed admirably. The method also 
proved useful when applied to a stratified flow/film cooled inlet nozzle boundary. 

The reference plane method of characteristics can be used on solid surfaces. While 
similar to the specification of normal pressure gradient under most conditions, the bound- 
ary condition appeared superior for the transient stability of the flowfield. 

Surface roughness can have an appreciable effect on the coolant requirements of a 
transpiration cooled nozzle. While the data from the test case examined were rather 
incomplete, it can be reasonably concluded that the design and fabrication of the transpira- 
tion coolant system employed coolant slots and channels that were far too coarse for the 
extremely thin boundary layers that exist on propulsive nozzles. Consequently, all future 
transpiration cooling systems need to consider the potential of significantly reduced per- 
formance due to surface roughness considerations. 

Coring of the centerline region of the nozzle can be used to restrict flowfield compu- 
tations to regions of physical significance. In the 1030:1 high area ratio nozzle analysis, 
the centerline region was only 20 % of the length of the wall and allowed significantly 
greater packing in regions of importance. 
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SECTION 9 

RECOMMENDATIONS 

As with the development of almost all good engineering tools, areas requiring further 
investigation increase with the greater understanding of the real flowfield requirements. 
Besides the perpetual need for greater documentation and demonstration cases, and the 
need to extend the modelling to three-dimensional systems, specific improvements of the 
Tethys model in the following areas are recommended and are discussed in greater detail: 

• Enhanced turbulence modelling 

• Inclusion of radiation heat transfer 

• Extension to other gaseous propulsion systems 

• Enhanced characteristic boundary conditions 

• Enhanced Neumann boundary conditions 

• Enhanced artificial diffusion modelling 

• Improved flowfield conservation 

• Reduced computational requirements 

Other areas of chemical propulsion consideration, such as solid particulate or liquid 
droplet interactions, are features that may be better suited to other computer models. 

Enhanced modelling of the fluid flowfield can be realized by the inclusion of more 
advanced turbulence models. The k-E equations for modelling turbulence are currently 
available in the newest version of Proteus and have been pasted into the Tethys code, 
awaiting implementation by the inclusion of appropriate calling, reading, and writing 
commands. Consequently, a more advanced modelling of the turbulent viscosity is readily 
obtained. However, the real arena in which turbulence modelling of a multi species 
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chemical rocket engine needs to be advanced is in the adequate modelling of the turbulent 
Prandtl and Lewis numbers. Throughout the demonstration studies, both the turbulent 
Prandtl number and the turbulent Lewis number were set equal to a constant, which was 
obtained by a ‘best guess’ type of engineering reasoning based on experimental measure- 
ments of a similar situation. More sophisticated models assume that the turbulence is a 
function of the laminar Prandtl number and the turbulent viscosity. However, in the pres- 
ence of more complex boundary conditions, such as surface roughness, the turbulent 
Prandtl number may be several orders of magnitude in error, since the application of the 
surface roughness model does not decompose drag into the components of skin friction 
and form drag, thereby invalidating the relations between momentum and heat transfer. 
Further, the turbulent Prandtl number may be significantly different in the transpiration 
cooled region of a rocket nozzle because the viscous transport modes of heat transfer are 
comprised of interdiffusional and Dufour heat transfer terms in addition to traditional heat 
transfer by conduction. In a similar manner, the turbulent Lewis number is undoubtedly a 
function of the strength of the Soret thermal diffusion term in addition to the gradients of 
all of the species present Since thermal diffusion may result in species gradients even 
along solid wall boundaries, the need for improved modelling of the turbulent Lewis num- 
ber is a more encompassing requirement. Consequently, it is the belief of the author that 
primary effort on turbulence modelling of chemical rocket engines focus on the unique 
aspects of heat and mass transfer in a chemical rocket engine while making use of the best 
available turbulent viscosity models. 

Inclusion of a radiation heat transfer model would enhance the predictive capabilities 
of both nozzle performance and more significantly, nozzle heat transfer. While this is pri- 
marily true for the combustion chamber region of the nozzle, some amount of radiative 
heating of the nozzle wall and other cold regions will occur throughout the converging- 
diverging portion of the nozzle and the fluid dynamics model should include this poten- 
tially important heat transfer mode. 

Currently, the Tethys model is capable of modelling the flow characteristics of nearly 
any ideal gas system with the caveat that, unless the fluid is comprised of a hydrogen/oxy- 
gen system, the fluid can not be subjected to any chemical reactions or species diffusion 
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and is limited to a maximum of six chemical species. Since this is a rather serious restric- 
tion, the model should be extended to at least eight chemically reacting and diffusing spe- 
cies equations. 

Enhanced characteristic boundary condition modelling requirements include 
improved methods of modelling flow conditions at fluid comers and along the centerline 
and improving the reasoning involved in the establishment of the boundary conditions 
needed. Improving the flow solution along the comer points of the flowfield offers the 
promise of increased numerical stability and accuracy. Enhanced centerline modelling is 
also critical since geometries that contain strong shock waves near the centerline are 
incompatible with the constant entropy assumption. The reasoning involving the switch 
from a characteristic relation to the specification of the static pressure at the outflow 
region of a nozzle also needs to be improved in order to more accurately assess the perfor- 
mance and heat transfer characteristics of a nozzle near the exhaust region. These aspects 
of improved boundary condition modelling will all be discussed in greater detail in the 
paragraphs to follow. 

Centerline boundary conditions need to be improved to enhance the ability of the 
Tethys model to predict shock waves and chemical reactions. Since the centerline of the 
nozzle represents a fluid pathline, a new set of characteristic relations needs to be devel- 
oped for application in this region, along with a means of implicitly formulating the spe- 
cies source terms along this pathline. 

Currently, the Tethys model employs all characteristic relations when the fluid contra- 
varient velocity (normalized) is greater than the speed of sound, but uses one Dirichlet 
boundary condition (specification of static pressure) when the contravarient velocity is 
less than the speed of sound. In reality, this is a distortion from what is actually appropri- 
ate since the static pressure should only be specified when the real fluid velocity is less 
than the speed of sound. In regions where the contravarient velocity is less than the speed 
of sound but the actual fluid velocity is not, a new boundary condition to replace that cur- 
rently being applied needs to be investigated. An alternative characteristics method, the 
bicharacteristic method, does not have this analytical shortcoming and should also be con- 
sidered. 
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Neumann boundary conditions, as applied to the solid or transpired surface of a 
rocket engine, should be extended to be conservative. Due to the complexity of a chemical 
rocket engine analysis, maximization of solution accuracy with a minimal amount of com- 
putational grid points is highly desirable. Changing to a conservative boundary condition 
formulation offers the promise of improved solution accuracy and reduced computational 
requirements. 

Artificial diffusion modelling in the Tethys code, which was simply extended from 
the Proteus code, was shown in the validation studies to be nonconservative for a radial 
flow system. Since these relations were undoubtedly developed for a cartesian coordinate 
system and then applied to an axisymmetric coordinate system, this finding is not surpris- 
ing, although it is very bothersome. Since the computational effort necessary for the calcu- 
lation of the artificial diffusion terms is essentially negligible, far greater emphasis 
directed towards complete inviscid flow conservation should be pursued. 

Further conservation properties of the flowfield should also be pursued. The current 
formulation of the flow variables is designed such that all of the gradient terms are conser- 
vative when the flow variables of concern are constant. In the authors terminology, this is 
only a first-order conservative system, whereas a second-order conservative system exists 
when the flowfield is completely satisfied for a linear property gradient existing in the 
flowfield. If this condition were satisfied, Couette flow modelling would be exact. That is, 
in the same exacting manner that the author insured that some shear stress terms were 
equal to zero in a Couette flow condition, the same exactness needs to be applied to the 
inviscid terms to insure that they are also equal to zero in a Couette flow simulation. 

The computational resources required of a complete Tethys flowfield simulation were 
extensive. Advanced computational techniques, such as multigrid acceleration techniques, 
and Runge-Kutta time step marching methods, should be considered as a means of 
improving the flowfield convergence rate. Likewise, reformulation of the implicit opera- 
tors, such as a linearized flux splitting method and/or decoupling of the species transfer 
equations from the continuity, momentum, and energy equations, may significantly reduce 
the computational requirements of a chemically reacting nozzle flowfield problem. 
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The methodology employed for the solution of axisymmetric flow problems in the 
Proteus code (Reference 5) is presented in this appendix, which is a reproduction of 
Appendix B in Reference 5. Although nomenclature and equation numbering represent 
those employed in the Proteus report, they are consistent with those used throughout the 
present work. The references cited within Appendix B of the Proteus manual are cited at 
the end of this appendix. It should be noted that the Proteus model has the capability of 
modelling swirling flow in the axisymmetric coordinate system, which is a feature not 
considered in the present work. 
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B. AXISYMMETRIC ANALYSIS 

The analysis used in Proteus for axisymmetric flow is essentially the same as for two-dimensional planar flow, 
described in the main body of this report. However, there are some additional terms in the axisymmetric equations 
that complicate things somewhat. For that reason, the axisymmetric analysis is described separately in this 
appendix. 

B.l GOVERNING EQUATIONS 

In cylindrical coordinates, the governing equations for axisymmetric flow, with swirl, can be written using 
vector notation as 


9(rQ) 9(rE) 9(rF) 9(rE v ) 9(rF v ) TI 
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Equation (B.l) thus represents, in order, the continuity, x-momentum, /■-momentum, ^-momentum (swirl), 
energy equations, with dependent variables p, pu, pv, pw, and E T . Note that the additional terms in these axis; 
metric equations destroy the strong conservation law form of the two-dimensional planar equations presented in J 
tion 2.1. Unfortunately, the axisymmetric form of the equations cannot be put into strong conservation law f< 
(Vinokur, 1974.) 


The shear stresses and heat fluxes are given by 
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dw 

r - =/i a7 



w 
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Rx 
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In these equations, x, r, and 8 represent the axial, radial, and circumferential directions, respectively; and u, v, 
and w represent the velocities in those directions. The remaining symbols are the same as those in the two- 
dimensional equations described in Section 2.1. 


For turbulent flow, /i, A, and k represent effective coefficients. The turbulence model is described in Section 
9.0. The only modification to the model for axisymmetric flow is the definition of IQI, the magnitude of the total 
vorticity. For axisymmetric flow. 
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When the generalized grid transformation of Section 2.3 (with y replaced by r), is applied to equation (B.l) 
the result is 
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Although this axisymmetric equation cannot be put into exact strong conservation law form, the procedure 
used to do so for the two-dimensional equation, described in Section 2.4, is nonetheless applied to equation (B.4). 
The result is 
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Using equations (B.2a) through (B.2g) these can be expanded as 
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B.2 LINEARIZATION 

Solving equation (B.5) for BQ/dr (assuming r is not a function of time) and substituting the result into the 
time differencing scheme of Beam and Warming, given by equation (3.1), for 3(AQ")/3r and 3Q"/dr yields 
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This equation must be linearized using the procedure described in Section 4.0. 
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B.2.1 Inviscid Terms 

For the inviscid terms the Jacobian coefficient matrix dEJdQ is 
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where f\ = u£, x + v£ r and / 2 = (£7- + p)lp. The Jacobian matrix 3F/3Q has the same form as 3E/3Q, but wit 
replaced by 7. 

For the additional term H, the linearization procedure gives 
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B.2.2 Viscous Terms 

To linearize the viscous terms, E Vi , E^, etc., must first be rewritten in terms of the dependent variables, 
with derivatives in the cylindrical coordinate directions transformed to derivatives in the computational direct 
using the chain rule. The shear stress and heat flux terms, given by equations (B.3) and (B.7), become 
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The above expressions for the shear stress and heat flux terms are substituted into equations (B.6e) through 
(B.6g). As in the two-dimensional planar case, the cross derivative terms are separated from the non-cross derivative 

I 

terms. In addition, for the axisymmetric case the non-derivative terms are included with the cross derivatives. 

The resulting five elements of E V[ (excluding the 1 IJRe, coefficient) are 
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+M$x(?rV«{ + Zx VV i + €i ww {) + VGAtrUUf + €x“V( + Zr^Wf) + ~ (£ x + f?)T{ (B. 1 le) 


For linearization it is convenient to rewrite the last element as 
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(Ev, ) 5 = + (M + A)£j{ir(uv) ( + k$ t y (£,. v 1 + < f x UV ) 


■ | [^ 2 (v 2 + w 2 )< + f r v + W 2 )«]+ J- (£ 2 + fi)T( 


The elements of Fy, have exactly the same form as those of E V| , but with £ replaced by tj. 
The five elements of Ev 2 (again excluding the 1 URe r coefficient) are 

(EV,)i =0 


(Ev,>2 = 2 HZ,T)xUn + 


1 


7x«n + -7 r(rv)„ 


+ u€r(TlrUn + r h v n) 


(Ev 2 ) 3 = 2 n£ r T\ r v n + k£ r 


1 


7x«n + - 7r(")n 


+ u^x(.n r u n + TJ X V„) 


W 


(Ev,)* = MtxVxWn + Mtr 7rW, - M$r - 


(Ev ,)5 = 2 M(.txlxU“n + Zr 7 rW,) + Af x 


1 


7x“«, + - nMrv) n 


+Ad 


7 x vu, + - 77 r v(rv) fl 


w 2 * 


(B.l 


(B.l 


(B.l 


(B.l 


(B.l 


+Mtx(nr VK„ + 7x VV , + 7 X W ,) + fl^irj.UU,, + > 7 X UV , + J7r tw,) -Mtr — + -j^r + ^rVr)T n (B. 1 


The last element can be rewritten as 

(Ev 2 )s = 2 M^xHxUU^ + 4,tJr w,) + ^x(7x«M, + + *4r<Jlx™ n + Hr w„) + A?, ^ (£ x u + 

w 2 /r 

+^x(7rV«, + 7x W, + TJ x WW n ) + ^4 r (T] r UU n + 77 X UV, + 7], WW q ) - fj£ r + — (£ x 77 x + ^ r T} r )T q (B.l 

r rr r 

The elements of have exactly the same form as those of Ev 2 , but with f replaced by 77 and tj replaced by |. 
The five elements of Hy are 


0 

II 

(B.l 

(Hv>2 = 0 

v kr 1 

(B.l 

(Hv) 3 =-2^ -- k(f z U' +rf t u„)+ -|^ r (rv). + 7 ,(rv),J 

(B.l 


169 


(Hy ) 4 = p(Zr W ( +T} r W n )-fJ — 


(H v )s = 0 


Performing the linearization, the Jacobian coefficient matrix 9Ey | /dQ is 


dEy, _ J_ 
3Q Re r 


r dE v ' 


dQ 


0 

5 (\ 


Si 1 


34 * 
r dE v , ' 


V 30 Al 
r dE v ^ 

V * A, 


df VP 

a n 


- a z 


0 

a m 


a. 


a n 


, , + a xr-r ( 

KPJ P 


- a 


dZ[pJ rr dZ[p 


+ a rr - r t 




' Sty, '' 

V * hi 


0 
0 

a (\ 


a r 


df [p 


' dEy ^ 


<*0 


a fdr' 

a^U e T/ 


where 


a xx = (2p + &)Zl + uZ 2 r 


a rr = tf 2 x+(2M + »Z 2 r 


a a =uZ 2 x+PZ 2 r 


a xr = (u + WZxZr 


a’xr = - ZxZr 

r 


' _ " cZ 
a rr ~ hr 

r 


<*o = 1 r(f 2 x+& 
Pr r 


f dE Vl ' 

dQ 
V / 2 I 


afa^j 3 f v l , v 

= ~a u — I - I- a„ — | - I -a„ — r ( 


dZ\PJ ” d$\p 


(B.13d) 

(B.13e) 


(B.14) 
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The Jacobian coefficient matrix for the remaining non-cross derivative viscous terms, 3F V| /3Q, has the st 
form as dE v JdQ, but with £ replaced by rj. 

And finally, linearizing Hy, the Jacobian coefficient matrix 3H V /3Q is 



where 


ORIGINAL PAQ€ * 
OF POOR QUALITY 
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B.23 Equation Of State 

The equation of state given in Section 4.3 must be modified slightly to add the swirl'velocity w. Thus, 


p = {y-\iE T 


■ ~ ^ p[u 2 + v 


2 +w 2 


or, in terms of temperature, 


T= — 


^- I |« 2 + v 2 + , 2 


P 2 


The derivatives arising from the linearization are the same as those presented in Section 4.3, except for 

dp y - 1 ( 


dp 2 


u 2 + v 2 + w 2 


dp 

d{pw) 


= -(y- D>v 


dT_ __ j_ 
dp c v 


£ T I I 2 2 2 

— I u 2 + V + w 2 

p 2 p 


(B.16) 

(B.17) 

(B.18a) 
(B.18b) 
(B. 18c) 
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dT _ w 
d(pw) c v p 


[f constant stagnation enthalpy can be assumed, the appropriate equation of state is 


r - 1 . 

p = — p\ 


h T - ^ (u 2 + v 2 + w 2 ) 


(B. 


and the temperature becomes 


T 



+ v 2 + tv 2 ) 


Again, the derivatives arising from the linearization are the same as in Section 4.3, except for 


dp 

dp 


r - 1 


h T + - ( u 2 + v 2 + w 2 ) 


dp y-\ 

- = tv 

d(ptv) y 

dT 1 / 2 2 2 \ 

— (ir + v* + tv *) 

dp c p p 

dT _ w 
d(p w) c p p 


B.2.4 Linearized Governing Equation 

The linearized form of equation (B.8) can now be written as 


(B. 


(B.2 


(B.2 


(B.2 


(B.2 



OWGINAL PAGE ft 
OF POOR QUALITY 



173 


AQ + 


0.A r 1 d 


l +&2 r 1 ^<9 


3E ’. .... 

n 35 |AQ 


a 

+ dr, 


fdt\" .. 

a *] 40 *feM 


0i Ar 1 
- < 


1 + 0? r 


_a_ 

df 




9E 


3Q 

Ar 1 f9(rE) 3(rF) 


AQ 


a 

+ d? 


/ . \" 
^3F, ' 


dQ 


AQ 


+ bg 1 40 


1 + 0 2 r l df dr, J 1 + 0 2 r 


At 1 


a(rE v ,) aCrF^) . 

df dr, 


V 


(I + 0 3 )Ar \_ 
l+0 7 r 


d(rE Vz ) dirty,) 
df + dr] 


Y* 


0 3 Ar l 
1 +0 2 r 


3(rE v ,) a(rFv 3 ) 

df + d n 


v 1 


AQ"'' + O ( 0 , - i - 0 2 j(Ar) 2 , (0 3 - 0,)(Ar) 2 , (Ar) 3 


(B.22) 


B.3 SOLUTION PROCEDURE 

Letting LHS(B.22) represent the left hand side of equation (B.22), we can write 


LHS(B.22) = 


K + - 


0]Ar 1 


1 + #2 r 


d_ 

df 


dE dE Vl 
3Q 3Q 


f 


dr/ 


dF dFy 
"dQ " dQ 


AQ" 


(B.23a) 


0,Ar 1 fd H 3H V 


where 

K = 1 + - 

l+0 2 r V^dQ 3Q 

and I represents the identity matrix. The term in braces in equation (B.23a) can be factored to give 


LHS(B.22) = 


K + 


0iA r 1 3 


l+0 2 r df 


dE dt Vi 
r — - - r 


3Q dQ 


I + K 


-i 


0, Ar 1 3 


1 + 0 2 r dr] 


dF 3F. 

r — = — r ■ 


/ Ar f 1 „-i 

a 

3E 3Ey 

a 

( \T 

3F 5F Vl 

li + Jr> K 

df 

'dQ r dQ J 

drj 

r — -r — +- 

^dQ dQ J 


dQ dQ 


AQ 


(B.23b) 


AQ 


(B.24) 


The last term represents the splitting error. 

Equation (B.22) can thus be rewritten in spatially factored form, and, neglecting the temporal truncation and 
splitting error terms, becomes 
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K + 


0 t Ar I 9 


( 


9E 96^! 

r 90 _r ~90" 


1 + 0 2 r 9£ 

Ar 1 f9(rE) 9(rF) 



K 


-i 


K + 


0>Ar 1 9 


1 + 0 2 r dq 


9F_ 9F V| 
r dQ~ r ~d^ 


AQ = 


I + & 2 r [ 9£ dq 


V 

+ — + — - + H + 

j 


Ar 1 


1 + 0 2 r 


9(rEy ) 9(rF v ,) . 


Y» 


(l+0 3 )Ar 1 
l+0 2 r 


3(rE V] ) 9(rF V; ) 
9<? + 9/7 


\" 


0 3 Ar 1 
1 + 0 2 r 


dq 


+ Hw 


9(rE v? ) 9(rF Vj ) 

+ ft; 


yt-l 


+ r-^~ AQ" _1 


l+0 2 


(B 


Using the procedure of Douglas and Gunn (1964), as written by Briley and McDonald (1977), equation (B.25) 
be split into the following two-sweep sequence. 


Sweep 1 (if direction) 


K + 


0i Ar 1 9 


1 + 02 r 9f 


9E 9E V , 

r — — - r 


9Q 9Q 


AQ =- 


Ar 1 ( 9(rE) 9(rF) 


1 + 02 r \ 9£ 9/7 


+ H 


Ar 1 
1 + 02 r 

0 3 A r 1 
1 + 0 2 r 


9(/-Ey t ) 9(rF y,, ) . 

+ + tlu 




9/7 


(1 +0 3 )Ar 1 

9(rE Vj ) 9(rF v , 5 )'| 

+ 1 + 0 2 r 

9^ ' ■ dh 

K J 


9(rE Vj ) 9(rFV 3 ) 

9f + ^ 


fl-1 


02 - n- 1 

+ — — AQ 
1 + 02 


(B.: 


Sweep 2 (/? direction) 


K + 


0)Ar 1 9 

1 + 0 2 r 9 q 


9F 

9Q 


— r 



AQ" = K"AQ* 


(B.: 


Or, expanding K and rearranging, 
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Sweep 1 (<* direction) 


-• 0 ( Ar 1 d 
AQ + 1 


AQ 


1 + #2 f df 

Ar 1 ^a(rE) a(rF) 
1 + 0 2 r ( d£ dr. 


0,Ar 1 d 
l + 0 2 r d£ 


d±, 


dQ 


v 


V'l. 

\ ) 


AQ* 


0,Ar 1 f3H 3H V Y -• 

+ - kx-^r UQ 


1 + 02 r 13Q dQ 


+ H 


Ar 1 


1 + fc r 


3(rE V| ) d(rF Vl ) 

df + dr, 


+ H V 


(1 + 0 3 )Ar [ 
l+0 7 r 


5(rE v ,) d(rF V: ) 

+ dr, 


0 3 A r 1 

1 +07 ~r 


ac^E^) + a(rF v ,) 

d$ + dr, 


V -1 


+ -^-AQ n ~ l 

1+07 


Sweep 2 (r, direction) 


AQ + 


0\Ax i a 
1 + 0 2 r dr. 



0| At 1 3 
l+0 2 r dr, 


f . v* 

ar 


dQ 


AQ 


0i Ar Ifd H 3HvT a* Ar 1 fdH 

1+02 r^3Q dQ) l+0 2 HdQ 



Applying the spatial differencing formulas of Section 5.0 results in 


Sweep 1 (<f direction) 
0,Ar 1 


AQ, + 


(1 + 0 2 )2Af r 


aftY A - ( aEY A * 

'^l, AQw Tad, 4Qw 


(1 + 0 2 )2(A^) 2 


~ [( r i-, fi-l + r ifi) K i1~\AQ,- 

-• 1 0, Ar 1 fdU dHv'] -• ' 

+2r i y} + r 1>1 / w r«?AQ i +(r,-/ i + r i+1 / >I )" i ? +l AQ i+1 J + j AQ, = 


Ar 1 


- S f (rE) + S n (rF) + H 


Ar 1 


l+0 2 

(1 +0 3 )Ar 1 ~ 
1 + 0 2 r . 


1+07 


. 1" 0 3 Ar 1 

^(rE^ + ^frF^J 


i[^(rE Vt ) + ^(rF V| ) + H v ]" 

±[* f (r£y,) + $,(/*,,,>]" ‘ 


AQ"' 1 


l+0 2 


(B.27a) 


(B.27b) 


(B.28a) 
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Sweep 2 (y direction) 


. * n 

AQj + 


0,Ar 1 


(l+0 2 )2A n r 


3F\ 

f 30i 


AQ 


7+i 


7+i 




0\ Ar 


(1 + 0 2 )2{Atj) 2 


rt%, 4Q"_, 


-<'h /<-. + *>/, + r M /„ + + 7 (If - ^ J*o; = 

l+«2 'Uq 3Q J 


(B.2 
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Appendix B 

Interior Point Numerical Algorithm 

A comparison between the solution methodology employed in the Proteus code and 
that employed in the Tethys code is described in this appendix. For the most part, the sim- 
ilarities between the models are extremely close, with the obvious exception that the 
chemically reacting and diffusing aspects of the Tethys model result in a significant 
increase in the fluid flow terms that need modelling. The other primary difference between 
the two codes is in the manner in which the grid transformation metrics and the radius are 
handled. A discussion of the Tethys code numerical procedure follows, with significant 
references to Appendix A (which is essentially Appendix B in the Proteus Analysis 
Guide). A discussion of the differences between the Proteus code and the Tethys code in 
the manner in which the mass, momentum, and energy equations are modelled is dis- 
cussed in Sections B.l and B.2. The equations of state are discussed in Section B.3, and 
the manner in which the species density equations are solved is discussed in Section B.4. 
The procedure used to specify fluid entropy and total enthalpy boundary conditions is dis- 
cussed in Sections B.5 and B.6, respectively. 

B.l. Finite Difference Approximation of Advective and Source Terms 

The spatial differencing in both the Proteus and Tethys codes employed central differ- 
ence approximations for both the implicit and explicit sides of the governing equations. 
Likewise, the grid transformation metrics were evaluated using central difference approx- 
imations in both the Tethys and Proteus codes, except at the boundaries where the Proteus 
code employed a three-point, second-order accurate approximation to the difference equa- 
tions normal to the surface while the Tethys code employed a two-point first-order accu- 
rate difference approximation. 

As discussed in Section 3, the governing equations used in the Tethys model were 
altered slightly from those used in the Proteus model by rewriting the radial pressure gra- 
dient into a weak conservation form, as illustrated in Table 20. Further, since the pressure 
term was rewritten in a conservative form, the K matrix presented in Equation B.23b of 
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Appendix A reduces to the identity matrix. 


Table 20. Comparison of the formulation of the radial pressure gradient. 


Proteus code 

Tethys code 

3(rP)-P 

8p 
r <5r 


Likewise, in the effort to maintain one-dimensional flow conservation, the evaluation of 
the radius in the two codes differed as presented in Table 21. 


Table 21. Comparison of the evaluation of the radius terms. 


Derivatives with respect 
to 

Proteus model 

Tethys model 

Radius 

Exact 

Exact 

Axial distance 

Exact 

Average 


In a similar fashion, the evaluation of the radius in the weak conservation form of the 
radial pressure gradient (in the Tethys code) was evaluated in an exact manner. 

B.2. Finite Difference Approximation of Diffusion Terms 
In an analogous effort to maintain one-dimensional flow conservation, the evaluation 
of the metrics and Jacobians in the viscous terms was modified from the methods used in 
the Proteus code. Due to the lengthy nature of the reevaluation of these terms, the evalua- 
tion of the viscous terms acting in the axial momentum equation will be elaborated on in 
detail, with the extension to the radial momentum equation, the energy equations, and the 
mass transfer diffusion equations being a rather straightforward extension of the concepts 
presented herein. 
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In the transformed coordinate system, the viscous terms of interest in the axial 
momentum equation are given in Equation (B.2.1) below (neglecting the reference Rey- 
nolds number): 

(B ' 2 ' 1) ^[j <^ xx + 4 r V] + ^[j <V*X + V„>] 

The normal and radial shear stress terms can be evaluated in a transformed coordinate sys- 
tem by applying the chain rule, as shown below. 

X 

( B .2.2) t xx = (2p + X) (^ x u^ + ri x u 11 ) + - [^(rvl^ + q^rv)^] 

An alternative formulation used in the Tethys model is shown below: 

(B ' 23) X xx = (2 ^ + ^ + W + + ^r V r| + 7 

For the purpose of comparing the finite difference procedures of the Proteus and Tethys 
codes, consider only the viscous terms that act on the axial velocity gradient: 

(B.2.4) (2(i + ^) (^ x u^ + r) x u T j) 

If the portion of the term given by Equation (B.2.4) is substituted into Equation (B.2. 1) 
and if only the £ difference portion of that equation is considered, the resulting expression 
is given by Equation (B.2.5): 

d f?ti i 7 ) ; d (2n i X t 

d£|_ J ( 2 ^ + ^) a ^J J ( 2 ^ + ^ aT1 _ 


(B.2.5) 
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In Equation (B.2.5), it is observed that the first term is comprised of ^-only deriva- 
tives acting on the axial velocity while the second term consists of ^ and T| derivatives act- 
ing on the axial velocity. The numerical approximation of Equation (B.2.5) is compared 
between the Proteus and Tethys codes on a dicretized finite difference basis with extensive 
reference to the arbitrary computational grid presented in Figure 83. In the four tables that 
follow, an examination of the procedures of handling both the noncross derivative shear 
stress terms and the cross derivative shear stress terms is presented. The comparison of the 
numerical approximation of other components of the diffusion terms is relatively straight- 
forward, with the fundamental differences between the two codes being in the manner in 
which the metric coefficients and the inverse Jacobians are calculated. 



Figure 83. Computational grid for comparison of viscous term finite difference 
approximations used in Proteus and Tethys codes. 
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The finite difference approximations of the £-only derivatives acting on the axial 
velocity of the axial shear stress term in the axial momentum equation are given in Table 
22 for the Proteus model. 


Table 22. Finite difference approximation of a £-only derivative component of the axial shear 
stress term in the axial momentum equation in the Proteus code. 
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The analogous approximation of the shear stress term in the Tethys code is presented 
in Table 23. 

Table 23. Finite difference approximation of a £-only derivative component of the axial shear 
stress term in the axial momentum equation in the Tethys code. 
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In Tables 24 and 25, the procedure for approximating a representative cross-derivative 
shear stress term is presented for the two different codes. 


Table 24. Finite difference approximation of cross-derivative components of the axial shear stress 
term in the axial momentum equation in the Proteus code. 
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Table 25. Finite difference approximation of cross-derivative components of the axial shear stres 
term in the axial momentum equation in the Tethys code. 


Tethys code approximation for 


d 

n 





It should also be noted that, when differenced near a boundary, the Proteus dicretiza- 
tion method described in Tables 22 and 24 replaces, where necessary, central difference 
methods with three point backward or forward difference methods for the evaluation of 
grid metrics and inverse Jacobians. 
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B.3. Equation of State 

Due to the multispecies nature of the flow model, coupled with the fact that a rocket 
engine is subjected to extreme variations in temperature throughout the flow field, the 
equations of state are significantly more complicated than those for the fully mixed, con- 
stant property, fluid considered in the Proteus model. As described in Section 2, the fluid 
temperature can be related to other fluid flow properties through the following relation: 

n 

P e - \ (u 2 + V 2 ) - p ^ C0j (« iR - C vi T R ) 

(B.3.1) -p _ i = 1 

Pl-Cvi 

i = 1 

and the pressure is related to the temperature and density through the ideal gas law: 


(B.3.2) 


P = pRT 


The relation of the fluid temperature to the dependent properties can be derived in the 
following manner: 

n / T \ 

(B.3.3) P e = X P i| W iR + J C vi dT 

> 


+ ^p(u 2 + v 2 ) 


i=l V t r 

Written in terms of only the fluid temperature and dependent properties. Equation (B.3.3) 
takes the following form: 


n-1 n- 1 T T j 2 

(B.3.4) Pe= XPi (M iR- W nR) + ZPij( C vi- C vn) dT + P w nR + Pj C vn dT+ 5-^- + 5-y 


i = 1 


i = 1 t r 


Integration of the Equation (B.3.4) and rearrangement yields the following: 

2 - . _ 2 n— 1 n- 1 

X Pi (W iR _M nR ) -p«nR+ £ Pj (C vi - C vn ) T R + pC vn T R 


2 p 2 p 


i= 1 


1 = 1 


n- 1 

^p.(C, i -C, n )+pC ¥n 


i = l 


(B.3.5) T = 
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While Equation (B.3.5) may give the appearance of being an explicit relation for the tem- 
perature as a function of the dependent variables, the temperature-averaged specific heats 
themselves, C v j and C vn , are a function of temperature, thus requiring an iterative solu- 
tion process in order to insure that the temperature is adequately satisfied by Equation 
(B.3.5). In the Tethys code, simple updating of the temperature averaged specific heats, 
following the solution of Equation (B.3.5) was used, and was able to converge to a solu- 
tion within 0.00001 % of the temperature within twenty iterations or less, usually being 
much less, and always requiring a negligible amount of computational effort. 

The derivatives of the temperature with respect to the other dependent variables can 
be determined by taking the derivative of Equation (B.3.4), as shown in Equation (B.3.6). 


(B.3.6) 


n- 1 n- 1 

d(P c ) = X (u iR"“nR) d Pi + X ( ^ vi_£vn) ( T ~ T R> a Pi + P C v aT 
i = 1 i = 1 


_ (u 2 + V 2 ) 

+ c vn (t - t r ) a P + « nR a P + u a (p U ) +v3(pv) — ap 

The dependence of the fluid temperature on the dependent properties is given by Equa- 
tions (B.3.7) through (B.3.11). 


(B.3.7) 

(B.3.8) 

(B.3.9) 

(B.3.10) 


9T _ (Cyi Cvn) (Tr T) ( M jR w nR^ 

apj pc; 

(T R -T)C vn ~u nR + i(u 2 + v 2 ) 

dp pc; 

5T u 

a(pu) " pc v 

8T v_ 

9(pv) “ pC v 

3T 1 

d(pe) pC v 


(B.3.11) 
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The relationship between the fluid static pressure and the dependent variables can be 
obtained by employing the relations presented above and realizing that the gas constant is 
also a variable. 


(B.3.12) 


P = pRT = 





M 

n 



T 


Using Equation (B.3.12), the dependence of the pressure on the dependent variables can 
be linked through the temperature as given by Equations (B.3.13) though (B.3.17). 


(B.3.13) 


3P_ 

3p. 


= RT 


L M i 


1 

M 


n dT 

+ pR 5— 

a P . 


(B.3.14) 


(B.3.15) 


(B.3.16) 


3P RT „ aT 

NT + pR ^ 


£ = 0 r*L_ 

a(pu) p d(pu) 


ap aj 

3(pv) P 3(pv) 


(B.3.17) 


ap = 3t 

3(pe) P a (pe) 


B.4. Modelling of the Species Continuity Equations 
The requirements to solve a set of species continuity equations in the Tethys model 
was performed by a straightforward extension of the procedure used to solve the continu- 
ity, momentum, and energy equations in the Proteus model, with complete coupling of the 
species continuity equations to the fluid flow equations, which for a six species system, 
required the simultaneous solution of a banded, tridiagonal, 9X9 matrix. In general, all of 
the relational dependences in the species continuity equations to the other dependent vari- 
ables are a trivial derivation, the exception being the relationship of the species source 
term to the dependent variables, which is presented in Appendix D. 
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B.5. Modelling of the Fluid Entropy 

At some computational boundaries, the relation between the fluid entropy and the 
dependent variables was required as the specification of fluid entropy was used in place of 
the total pressure specification. The equation for the entropy of a mixture of species is 
derived in Reference 46 and is shown in Equation (B.5.1) below: 

n n 

(B.5.1) S = R - Rln ( p ) - R X x i ln(x i } 

i = 1 1 i = 1 


where the nondimensional entropy (i.e., s/R) of species i is determined from a polynomial 
curve fit (Reference 47) and is based on conditions at 1 atmosphere. Likewise, in Equation 
(B.5.1). the units on the pressure are in atmospheres. 

In all of the test cases considered, the boundary conditions on the species mass frac- 
tion were similar to the those used on the entropy. That is, in cases where the Dirichlet 
boundary conditions were used on the entropy, Dirichlet boundary conditions were also 
used on the species boundary conditions. Consequently, the gas constant will follow the 
variations in the fluid entropy and the following operation can be performed to determine 
the relation between the entropy and the fluid temperature and pressure. For Dirichlet type 
boundary conditions on the entropy, the following expression can be written: 


S 

= <S> 




At= (£) 


i 3 t i ap 


where the relation between the pressure and the temperature to the dependent variables is 
given in section B.3. 

For cases where the computational gradient in entropy and species gradient is equal to 
zero, the following relations can be used: 


(B.5.3) 


3s 


R3(^> 


+ =0 

R&n 


where an analogous relation exists if applied along a £, line. Making use of the fact that the 
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gas constant is a function only of the mass fractions of each species. Equation (B.5.3) sim- 
plifies to the following boundary condition that was modelled in the Tethys code: 


(B.5.4) 


3 <s _ 

dn 


Equation (B.5.4) can be written in the following implicit manner: 


(B.5.5) 
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where the relation of the temperature and pressure to the dependent variables is given in 
Section B.3. 


B.6. Modelling of the Total Enthalpy 

The specification of the fluid total enthalpy, in either Dirichlet of Neumann form, was 
often employed in the modelling of the boundary conditions along a rocket nozzle, with 
inlets having a fixed (i.e., Dirichlet) boundary condition while the gradient in total 
enthalpy (i.e., a Neumann boundary condition) was often employed along a slip wall or a 
centerline. The equation for the total enthalpy is given in Equation (B.6.1): 


(B.6.1) 


pe + P 

P 


For Dirichlet type boundary conditions where the total enthalpy desired is specified, 
the following implicit relations are maintained: 


(B.6.2) 




At 


where the time change of the total enthalpy can be derived from Equation (B.6.1) as 
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written as shown in Equation (B.6.3): 


(B.6.3) 
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where the relation between the static pressure, P, and the dependent variables is given in 
Section B.3. 

In cases where Neumann boundary conditions are employed, the following relation- 
ship is employed: 


(B.6.4) 
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pe + P 
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where the computational gradient given by Equation (B.6.4) equals zero. An analogous 
expression arises if one specifies the gradient in total enthalpy in the £ direction. 
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Appendix C 

Multicomponent Species Diffusion 


The laminar diffusion of a multispecies gas can be contributed to a gradient in species 
concentration or a gradient in temperature or pressure. The contributions can be written as: 


(C.0.1) 


j, = j , 00 +j, rn +j< p) 


where jP^is the diffusion flux as the result of species concentration gradients, is the 
diffusion flux due to temperature gradients, and is the diffusion flux due to a pressure 
gradient. The components of diffusion have the characteristic that the sum of the thermal 
diffusion species flux is equal to zero. Therefore, since diffusion as the result of a pressure 
gradient is assumed to be very small in a propulsive nozzle, the following equations for 
the components of species flux must be true: 


(C.0.2) 


n 



i = 1 


V i< T > - n 

(C.0.3) 2 j 1 i " 0 

i = 1 


In the sections to follow, the evaluation of the individual components of mass diffu- 
sion will be discussed. 


C. 1 . Multispecies Ordinary Diffusion 

As described by Bird et al. (Reference 8), the multispecies transfer equations of ordi- 
nary diffusion can be written as: 



2 

P 1 J !J 1 

j = l 


(C.1.1) 
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In general, the multicomponent diffusion equation given in Equation (C. 1 . 1 ) is rarely used 
directly due to the complexity of the multicomponent diffusion term, D-j, which is, in gen- 
eral, a function of both temperature and species concentration. Thus, experimental evalua- 
tion of such a term would require a large number of variables and is consequently rarely 
pursued. 

Instead, the Stefan-Maxwell equations of multicomponent species transfer, which 
relate multicomponent species transfer as a function of binary diffusion coefficients, 
are the usual starting point for ordinary, multispecies, mass transfer calculations. 




p ^ 

j = i 
j*i 


j= 1 J 


Rigorously, the Stefan-Maxwell equations provide a relation for the multicomponent 
transfer of species in a first-order sense only. However, extensive experimental tests and 
analytical predictions have shown that second-order effects, as concerns ordinary diffu- 
sion, are quite small (Reference 48). 

The significant advantage of the Stefan-Maxwell equations is that they algebraically 
relate the multicomponent ordinary diffusion fluxes to the binary diffusion coefficients. To 
a first-order approximation, the binary diffusion coefficients are independent of species 
concentration and are amenable to experimental determination. 

Some simplification in the calculation of the diffusion fluxes can be made by making 
use of the fact that the sum of the ordinary diffusion fluxes is equal to zero. 
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Further simplification to the right-hand side of the above equation could be made through 

n- 1 

the use of the mass conservations relation, co n = 1 - ^ C0j . However, the determination 

j = l 

of the multicomponent thermal diffusion coefficients requires the knowledge of the rela- 
tion between ordinary species diffusion and all of the species concentration gradient 
terms. For a four species system, Equation (C.1.3) can be written in the following matrix 
form: 



^11 4*12 ^13 


j,® 


0.. 


04. 

(C.1.4) 

^21 ^22 ^23 
_ ^31 4*32 ^33. 


.,x, 


P.2 

P» 

VcOj + ... + 

042 

043 


The extension to other species systems is apparent Equation (C. 1 .4) is solved by fac- 
toring the left-hand side into lower-upper triangular matrices and then successively solv- 
ing for the species dependence vectors, (5. . The resulting solution is in coefficient form: 

(C.1.5) ji X) = i C ij Vt0 j 

j = l 


Employing the mass conservation rule. Equation (C.1.5) can be put into a relation 
involving n - 1 species: 


(C.1.6) 



n- 1 


XV 

j = l 


( 0 . 

J 


where A- = - C in . 


C.2. Multispecies Thermal Diffusion 

The determination of multispecies thermal diffusion is neither a well-developed nor a 
simple methodology. In contrast to the calculation of multicomponent ordinary diffusion 
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where the Stefan-Maxwell equations allowed the determination of multicomponent diffu- 
sion based on binary diffusion coefficients, the complex nature and relatively poor devel- 
opment of the multicomponent thermal diffusion factor forces reliance on purely 
theoretical calculations. Given these conditions, an expression developed by Waldman 
(Reference 49) was adopted for the calculation of the thermal diffusion coefficient. Wald- 
man gives the following expression for the multicomponent thermal diffusion factor: 


(C.2.1) 
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a ij = 31M 


i + MjJ L2 M. 


The coefficient can be determined using the following first-order approximation: 
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(C.2.2) 


Z x j (a ij A i +a ii A j ) - T 

j=l 


The ajj and a- j coefficients are a function only of the binary masses and the temperature 
of the species. 

M i M j 

(C2. 3 ) a ij = 5 M i + Mj F ij 


(C.2.4) 
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where the coefficients F- • and F-j are defined as follows: 
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(C.2.6) 
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The nondimen sional terms. A— , Bjj, and C-, are functions of the Q integrals, which in 
turn are functions of the properties of the species-pair combination and the fluid tempera- 
ture (Reference 48). 

Using the expressions presented in Equations (C.2.3) to (C.2.6), Equation (C.2.2) can 
be solved by the simultaneous solution of a set of equations equal in number to the number 
of species present. Once Equation (C.2.2) has been solved. Equation (C.2.1) can be solved 
for the multicomponent species thermal diffusion factor, ou . As defined, the thermal dif- 
fusion coefficient, Dj , is related to the diffusion factor by the following relation: 


(C.2.7) 


D J = Z°U (k T>j = X D ij X i X X k“j 

j = l j = l 


T 

jk 


k = 1 


where Djj is defined as the coefficient of ordinary species diffusion when based on species 
mole fraction: 

(C.2.8) 


j = i 


A relation for the multicomponent diffusion coefficient based on mole fraction, Dy, 
and one based on mass fraction, C-, can be determined by equating the flux relations: 


(C.2.9) 




( 0 - 

J 


Substitution of the relation between species mass fraction and species mole fraction can be 

* 

made in Equation (C.2.9) to obtain the relation between Djj and Cjj. 
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T 

The final form of the equation for the thermal diffusion factor, Dj , shown in Equa- 
tion (C.2.7), was used for the evaluation of the Soret mass and Dufour heat transfer coeffi- 
cients. Equation (C.2.10) was used for the evaluation of the ordinary diffusion coefficient 
required in Equation (C.2.7). 
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A ppendix D 

Species Source Function 


The species source function, a., is a function of species concentration, reaction 
mechanism, and reaction rates. In Reference 50, the species source function is derived as 
presented below: 


(D.0.1) 


< 3 . 

l 


m 

= M. V Av.. 
i Lu lj 

j = l 
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V.. 







where the terms are defined as follows: 


rr 

v y = Product stoichiometric coefficient of species i for the jth reaction mechanism 


v.j - Reaction stoichiometric coefficient of species i for the jth reaction mechanism 


Av.. = v.. - v.. 
ij U iJ 


K f j = Forward reaction rate for the jth reaction 


K b j = Backward reaction rate for the jth reaction 


The reaction mechanism and the forward reaction rates for the gas system of interest 
(hydrogen-oxygen) are presented in Appendix I. 

The forward and backward reaction rates can be related to the equilibrium constant, 
Kpj, through the following expression (presented in Reference 23, p. 22): 


(D.0.2) 


Kr: _ -Av. 

^ = K pj (RT) J 
K bj J 


n 

where Av. = Y Av.. 
J ^ ij 
i = 1 


is the net mole production of the jth reaction. The equilibrium 
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constant can be determined as a function of the stoichiometric coefficients and the Gibbs 
Function, g-: 


(D.0.3) 


K p j = exp 


-5>i A y (RT) 

i = 1 




J 


where the Gibbs function used in Equation (D.0.3) is expressed on a molar basis. 

Substitution of Equation (D.0.2) and (D.0.3) into Equation (D.0.1) yields the follow- 
ing relation for the species source function: 


(D.0.4) o. = 

j= 1 



The dependence of the species source function on the dependent variables can be 
determined by partial differentiation of Equation (D.0.4) with respect to the dependent 
properties and with the realization that the backward reaction rate, K b j, can be a strong 
function of temperature (Appendix I), and for catalytic reactions, K bj is also a function of 
the species concentration. 

D.l. Relation of the Species Source Function to the Dependent Variables 
To accurately describe the dependence of the species source function on the fluid den- 
sity and the species densities, it is necessary to recast the species source function of Equa- 
tion (D.0.4) into a form that is independent of the final species density term, p^, as this 
term is not a dependent variable since it is related to the other species densities by the rela- 
tion: 

n- 1 

(DM) Pn = f-Z p i 

i = 1 

Using Equation (D.1.1), terms contained in Equation (D.0.4) can be recast in terms of 
the dependent species densities solved in the Tethys model. Thus, 
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The derivative of Equation (D.1.2) with respect to the fluid and species densities is as fol- 
lows: 
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In a similar manner, the derivative of Equation (D.1.2) with respect to the species densities 
is as follows: 


(D.1.4) 




The complete dependence of the species source function with respect to the fluid den- 
sity is shown in Equation (D.1.5): 


do. 

roi ' 5 ) 


Av.. 
‘Z_r >J 

j = l 




i = 1 


K, 


bj 


+ M 


■ Y Av.. 

'jLd >J 


j = l 




MJ li'M- 
i = 1 ' i = 1 


<* bj 

5P 


+ Av.. 
'z-d >J 
j = l 


[- 4 v,K pj 1ST] * (ST)‘"® P ']n g) ” 

i = 1 _ 


K 97 



200 

The dependence of the species source function on the species densities has a similar form. 


da. 

(D.i.6) 
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The dependence of the species source function on the other dependent variables is 
significantly simpler, being a function only of the relation between the particular depen- 
dent property and the fluid temperature. 
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The dependence of the fluid temperature on the dependent properties is described in 
Appendix B, while the dependence of the backward reaction rate, K b j, and the equilibrium 
constant, K p j, are described in the following sections. 


D.2. Backward Reaction Rate, K b j 

The backward reaction rate of the hydrogen-oxygen system of interest is assumed to 
obey a modified form of the Arrhenius equation (Reference 23): 


(D.2.1) 


K bj = B/AjxpI-Ej/fRTH 


where Aj, otj, and Ej are experimentally estimated coefficients. When third body (i.e., cata- 
lytic) chemical reactions are considered, the assumption is made (supported by theoretical 
and experimental measurements, as discussed in Reference 1) that the reaction rate coeffi- 
cients aj and Ej are independent of the third body reaction. Consequently, a mean third 
body reaction rate, Aj, can be used, which is defined by the following expression: 


(D.2.2) 


K h = ^ 
bj 
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exp [-Ej/ (RT) ] 


where the mean reaction rate is equal to the density weighted sum of the third body reac- 
tion rate coefficients: 


(D.2.3) 
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To facilitate a simple derivation. Equation (D.2.3) can be recast in terms of the fluid 
density and the dependent species densities: 

n - 1 


a ; 


(D.2.4) K bj = T J exp [-Ej/ (RT) ] 
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In this form, differentiation with respect to the dependent variables is readily accom- 
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plished. The results are presented in Equations (D.2.5) to (D.2.11). 

For catalytic reactions, the dependence of the backward chemical reaction rate is 
shown in Equation (D.2.5) and (D.2.6). 


(D.2.5) 
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For binary exchange reactions, the backward chemical reaction rate is related to the 
species densities only in there relation to the fluid temperature, as shown in Equation 
(D.2.7) and (D.2.8): 
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For either catalytic or binary exchange reactions, the backward chemical reaction rate 
is independent of pu, pv, and pe, except for the dependence of these properties on the fluid 
temperature: 
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(D.2.11) 
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The dependence of the fluid temperature on the dependent variables is presented in 
Appendix B. 

D.3. Equilibrium Constant, K p j 

The relation for the equilibrium constant, in terms of the Gibbs free energy and the 
stoichiometric coefficients, is given by Equation (D.0.3), which is repeated here. 

n 

(D.3.1) Kpj = expf-^gjAVyAR/T) 1 

1= 1 ' 

where the Gibbs Function, gj , is defined as: 

(D.3.2) gi = hi-Ts? 


Here, the superscript o signifies evaluation at one atmosphere of pressure. 

The following form of curve fit to the species enthalpy and entropy is adopted from 
Reference 47: 
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Substitution of Equations (D.3.3) and (D.3.4) into Equation (D.3.2) and rearranging 
terms yields the following expression: 
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The dependence of the equilibrium constant on the temperature can be determined by tak- 
ing the natural logarithm of Equation (D.3.5) and differentiating: 


(D.3.6) 
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Further rearrangement and resubstitution yields the following expression for the relation- 
ship between the equilibrium constant and the fluid temperature: 


(D.3.7) 
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The equilibrium constant is seen to be a function only of the fluid temperature. The 
relationship between the fluid temperature and the dependent variables is given in Appen- 
dix B. 


205 


Appendix E 

Time-Averaged. Density- Weighted. Turbulent Flow Equations 

The use of a density-weighted technique has come into favor as the preferred means 
of time-averaging the Navier-Stokes equations in compressible flow situations. The aver- 


aging process is performed by first decomposing the velocity, species densities, tempera- 
ture, energy, enthalpy, and total enthalpy into mean and fluctuating components: 

(E.1) 

v = v + v" 

(E.2) 

w = w + w" 

(E.3) 

u = u + u" 

(E.4) 

T = T + T" 

(E.5) 

Pi = Pi + Pi" 

(E.6) 

e = e + e" 

(E.7) 

h = h + h" 

(E.8) 

H = H + H" 

where the tilde (~) above the properties is defined as a time-averaged, density weighted 
term divided by the time-averaged density. The tilde signifies that the following operation 
has been performed: 
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(E.9) 


X = 


pX 

7 


Since X = X + X" and pX = pA. + pA." , the averaging process yields: 


(E. 10) 


pA. = pA.+ pX" 


Division by the time-averaged density yields: 


(EM) 


pX 

7 




Since the time-averaged density, p, is a finite number, it must be true that the numer- 
ator of the second term on the right hand side of Equation (E.ll), , is equal to zero. 

This correlation is fundamental to the density-weighted technique that follows. It will be 
further assumed that the time-average of all diffusion terms is equal to the density 
weighted, time-average of each of the individual terms comprising the constitutive equa- 
tion. For instance, the time-averaged shear stress, f , is assumed to be equal to the fol- 
lowing: 


(E.12) 



The fluid viscosity, and all other properties, including chemical reaction rates and 
equilibrium constants, are evaluated at the time-averaged temperature and species concen- 
tration. 

Recall the cylindrical coordinate system form of the Navier-Stokes equations: 
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where the centrifugal force, pw 2 , is retained, even for two-dimensional flow with no aver- 
age tangential flow. For brevity of approach, the alternative to the description of the radial 
pressure gradient is shown in Equation (E. 13) and, as will be shown in the following dis- 
cussion, the turbulent aspect of the flow has no effect on the pressure gradient terms other 
than to change the terms to a time-averaged quantity. 

Decomposing the enthalpy( H = E T +P = p[h+ (u +v )/2] ), velocity, mass 
density, and energy into fluctuating and time-averaged quantities. Equation (E.13) can be 
written as shown in the following equation: 
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p(u + u") 
p(v + v") 
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(E.14) 
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terms 


Note that the density has not been decomposed into average and fluctuating components. 
This is an attribute of the density weighting technique. If the above equations are time- 
averaged, the density-weighted, time-averaged, Navier-Stokes, equations are: 
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The time- averaged terms containing the stagnation enthalpy, H", should be decom- 
posed into a more fundamental form: 


(E.16) H = p [h + (u 2 + v 2 + w 2 ) /2] 

Furthermore, the enthalpy, h, is the sum of the mass weighted species enthalpies: 

n 

(E.17) ph = X p i h * 

i = l 

Now, decomposing the stagnation enthalpy into mean and fluctuating components 
yields: 
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n 

(E.18) H + H"= ^ (pj + p/') (hi + h/') +Pt(u + u") 2 + (v + v ") 2 + w ,, 2 ]/2 
i = 1 

Multiplication of Equation (E.18) by the fluctuating axial velocity, u", and time-averaging 
results in the following expression: 

n 

(E.19) lTiP= (hipjV' + PihjV') + vpiTV' + upu" 2 + (pu" 3 + pv" 2 u" + pu"w" 2 ) /2 
i = 1 

Similarly, multiplying Equation (E.18) by the fluctuating radial velocity, v", and time- 
averaging yields the following: 

n 

(E.20) HV* = ^ (hiP^V' + PjE/V') + upiTV' + vpv" 2 + (pv" 3 + pv"u" 2 + pv"w" 2 ) /2 
i = 1 


Assuming that the triple correlation turbulent fluctuations are negligible, the stagnation 
enthalpy expressions in Equations (E.19) and (E.20) simplify to the following: 

n 

(E.21) IFu 7 ' = ^ (hjp/V'-f- pjh/'u") + vpu"v" + Gpu" 2 

i = 1 

n 

(E.22) !rV' = ^ (hipfV' + Pi^V') + 5fnPV' + vpv" 2 

i = 1 

In general, the time-averaged pressure, P, and species source function, cr, will be 
functions of some time-averaged turbulent terms (e.g., P = pRT + Turbulent terms). 
However, for the flow problems of interest in this study, these terms are small compared to 
their time-averaged counterparts and consequently were not considered. 
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Summarizing, the dominant density-weighted, time-averaged, Navier-Stokes equa- 
tions are: 




pu 

pu 2 + pu" 2 + P 
puv + pu"v" 

(h i p/ 7 u 7 ' + p^ r a') + vpuV + upu" 5 

Pi“ + Pj" U * 



Viscous 

terms 



211 


The equations of state are taken to be unchanged by the presence of turbulence: 


(E.24) 


P = pRT 


where the gas constant of the mixture, R, is equal to the universal gas constant divided by 
a mean molecular weight: 


(E.24.1) 


R = = 

M 


where the mean molecular weight of the mixture is defined as: 


(E.24.2) 


I 


The relation of the total energy to the internal energy and the velocity is also assumed 
to be unchanged by the presence of turbulence: 


n (T \ 

pe = ^pj JC v .dT + u iR +p(G 2 + v 2 )/2 
i = 1 Ur > 
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Turbulence Modelling 
F.l. The Boussinesq Approximation 

The time-averaging of the Navier-Stokes equations results in the creation of numer- 
ous turbulence quantities for which no known solution exists. These terms are, however, 
usually modelled by assuming they have ‘viscous-like’ behavior, commonly referred to as 
the Boussinesq approximation (Reference 51). For example, the turbulence term, p u 'V\ 
is usually modelled as follows: 


(F.1.1) 



In a similar fashion, turbulent heat transfer and mass fluxes are approximated by the 
following diffusion equations: 


(F.1.2) 



i = pp!m 


3o> i 

8r 


(F. 1.3) 


(F.1.4) 


P< h "v" = - P CpT"v" = k> ^ 

i = 1 


n 





v' - - , OUJ i 

= 

i = 1 


Analogous expressions exist for the turbulent stress terms containing the axial velocity. 

Turbulent shear stress terms containing the square of the velocity are also assumed to 
have viscous-like behavior analogous to the laminar counterpart: 
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Similar expressions occur for the modelling of pv"v" and p w"w" . 

Further analogies between turbulent stresses and their laminar counterparts can be 
made by introduction of turbulent Prandtl and Lewis Numbers. 


(F.1.6) 


Pr‘ 



(F.1.7) 


t P C p ^iM 
Le = - 

k l 


The expressions given by Equations (F. 1.1) to (F.1.7) are the typical starting points 
for turbulence modelling, with the modelling goal being to find suitable correlations for 
the turbulent coefficients introduced by the Boussinesq approximation (i.e., (I 1 , k l , , 
and X 1 ). Primary focus has been on finding a suitable correlation for the turbulent viscosity 
and modelling the remaining turbulence terms through either the introduction of a con- 
stant or an algebraic expression. In this study, the turbulent viscosity was modelled by the 
use of an algebraic turbulence model, whereas the turbulent Prandtl and Lewis number 
models were assumed to be constant. In the sections to follow, a description of the various 
models applied for this study or available to use in the Tethys flow modelling computer 
code, is made. 


F.2. Algebraic Turbulence Model 

The basic algebraic turbulence model adopted for this study is the Baldwin-Lomax 
model (Reference 52), which is a two-layer algebraic turbulence model. Thus, 


(F.2.1) 




^inner ^ or ^n~^b 

^ > 
u l „ for y„ > y, 
'outer J n J b 

w - 


where y n is the normal distance from the wall and y^ is the distance at which the outer 
and inner viscosity values match. 
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The inner region turbulent viscosity in the modified Baldwin-Lomax model has the 
following form: 


(F.2.2) 


I if = p/ 2 |Q| 

'inner r 1 1 


where / is the mixing length and Q is the total vorticity. The mixing length is given by the 
expression: 


(¥.23) 


l = Ky n [ 1 - exp (-y + / A + ) ] 


and the magnitude of the total vorticity is defined as: 


(F.2.4) 


|Q| = 


9v _ 3u 
3x dy 


In Equation (F.2.3) for the mixing length, the variable k is the von Karman constant 
(taken to be equal to 0.4), the variable y + is the nondimensional distance from the wall. 


and A + is the Van Driest damping factor: 
(F.2.5) 


+ w 

y = ~ir ~ y « 


(F.2.6) 


a+ 26 
A = — 2 

N 2 


Where N 2 is given by the following relations: 


(¥. 2 . 1 ) 


n2 = 1 +i1 4Jpam 3 


W 


3<t> 

Solid wall regions 
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(F.2.8) N 2 = - 


3P 


1 


(PV) T 


1 - exp 


,11.8 (pV) a 


w w 


w r w 




X 

w w 


+ exp 


,11.8 (pV) u 


w r w 




X 

w w 


Transpired wall regions 


where 3P/5<)> is the pressure gradient tangential to the wall. In the original Proteus 
2 

code, N was equal to unity. The modified form is attributed to Cebeci (Reference 53). 
In the outer region, the turbulent viscosity is computed from the following relation: 


(E2 * 9) ^outer _ ^ C cp F Kleb F wake 

where 91 is the Clauser constant, taken as 0.0168, and C r _ is also a constant, taken to be 
1.6. The wake parameter, F wake , is specified by the expression given in Equation 
(F.2.10): 

(F- 2 - 10 ) F wake = >Wx F MAX 

where ^max * s maximum value of the following correlation: 

(F.2.11) Fmax = y n |0|[l -exp(-y + /A + )‘ 

Wall bounded flows 

(F-2-12) F max = y n |£i| 

Free turbulent flows 

where y MAX corresponds to the location at which F MAX occurs, as measured from the 
wall for wall bounded flows and as measured from the location of the maximum or mini- 
mum velocity for free turbulent flows. The Klebanoff intermittency factor, F^t,, accounts 
for the intermittent turbulent nature of the flow at regions significantly far away from the 
wall. This factor is given by the following relation: 
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(F.2.13) 


F = 
r Klcb 


1+B 


^Klebyn ^j 6 

^MAX ) . 


where B and are constants equal to 5.5 and 0.3, respectively. 

A special exception to the formulation of the outer region turbulence model was 
employed for the analysis of a film cooled rocket nozzle, as discussed in Section 7.2. This 
change was performed so that a reasonable solution for cases having both wall bounded 
turbulent flow and free stream turbulent flow could be obtained. 


F.3. Turbulence Modelling for Roughened Surfaces 
The alteration of the turbulence model in consideration of a roughened surface was 
performed by using a ‘displaced distance’ model, as described by Cebeci and Smith (Ref- 
erence 54). In this model, the equations that relate the turbulent viscosity to the distance 
from the wall is modified by an expression that is a function of the surface roughness, as 
given below: 

(F.3.1) Ay + = 0.9^-^k^- k + exp (-k + /6.) J 

The nondimensional surface roughness parameter, k + , is a function of the surface rough- 
ness (in equivalent sand grain roughness measures), the shear stress, and the fluid density 
and viscosity: 

+ kp 

(F.3.2) k = 



In the application of Equations (F.3.1) and (F.3.2), since the wall is displaced by a sur- 
face roughness parameter, the turbulent viscosity is nonzero even at the wall. Hence, the 
wall shear stress as used in Equation (F.3.2) is a function of both the laminar and turbulent 
viscosities at the wall. 
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The effect of surface roughness on mass and heat transfer is usually not as pro- 
nounced as the effect on skin friction, as discussed in Reference 14. However, no correc- 
tion to the turbulent Prandtl number, Lewis number, or the boundary conditions to account 
for this phenomena was applied in the turbulence model. 


F.4. Turbulence Modelling of Pr\ Le\ and X 1 
The turbulent Prandtl Number was assumed to be either a constant, which may vary 
depending on the specific application, or to have the following relationship with the turbu- 
lent viscosity and the laminar Prandtl number: 


(F.4.1) 


Pr 1 


'Pr3 


CprjPr 


1 - exp 


' C Pr4 ^ 

■, pVpJ 


( 

1 - exp ■ 

V 


C Pr2 ' 

(Pr) pVpJ 


where the constants Cp^, C pr2 , C pr3 , C pr4 , are 0.21, 5.25, 0.20, and 5.0, respectively. 
Equation (F.4.1) was excerpted from the original Proteus code and was never applied in 
any Tethys fluid flow simulations. 

The turbulent Lewis number, Le l , was assumed to be constant. Due to the analogy 
between the transport of species and heat, the Lewis number is likely close to unity. Thus, 


(F.4.2) 


Le l = 1.0 


The turbulent second coefficient of viscosity, A. 1 , was assumed to be related to the first 
coefficient of turbulent viscosity, |i l , in the same manner as their laminar counterparts: 

^ _ A, _ _2 
M- 3 


(F.4.3) 
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A ppendix G 

Transformation of the Governing Equations 

The axisymmetric version of the Navier-Stokes equations modelled in this study were 
transformed into a generalized coordinate system. Transformation allows significant ease 
in modelling the flow equations by the use of central difference approximations, more than 
offsetting the complexity involved in the transformation. The general equations governing 
the flow of a multispecies fluid have the following form: 

(0.1) s («Q> + K (rE)+ 3? <*> + H - K <rE v ) + | (rF v ) + H v 

Consider the generalized transformation of the equations into the following general- 
ized coordinate system: 

(G.2) T = T(t) 

(G.3) % = % (x,r) 

(G-4) T| = T| (x,r ) 

Using the above transformations of physical space into computational space and assigning 
the time derivative to be unaltered by the transformation (i.e. x = t), the following chain 
rules can be applied: 


(G.5) 

90 

9!; 90 9r| 90 

= 

9x 9^ + 9x 9rj 

(G.6) 

3<t> 

9^ 90 9q 90 

Tr = 

9r 9^ + 9r 9r| 


Applying the generalized coordinate transformation from physical to generalized 
space, the governing flow equations yield the following result: 
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3 (rQ) + H k (rE " rE v) + ^ 4 (rE - rE J 


dz 


9x 8q 


(G.7) 


94 9 . . 9r| 9 . . 

+ 57 5l' rF - rF v) + g; ^(rF-rF v ) + H-H v = 0 


where the nonconservative pressure gradient term in the radial momentum equation takes 
the following form: 


(G.8) 


9P _ [94 9P 9q 9P~ 

r dr ~ T |_9r 94 + 9r 9q 


In the present form, the generalized coordinate transformation is in a nonconservative 
form since the transformation metrics are outside of the partial derivative operators. To the 
extent possible, the transformation should be put into conservative form, and the transfor- 
mation metrics should be expressed in terms of physically determined spatial coordinate 
differences. These operations are performed in the expressions that follow. 

Computational distances can be written in terms of the metric derivatives and physi- 
cal distances, as outlined below: 

(G.9) d4 = 4 x dx + 4 r dr 

(G.10) dq = q x dx + q r dr 

Similarly, physical distances can be written as a function of computational distances: 
(G.ll) dx = x^d4 + x^dq 

(G. 12) dr = r^d4 + dq 

Solving Equations (G.9) to (G.12) for d4 and dq yields: 


[ Vil" r ^ x q ]d ^ =r q tix " x q cir 


(G.13) 
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(G.14) [x^ - r^xj:] dr| = r^dx - x^dr 

Comparing Equations (G.13) and (G.14) with Equations (G.9) and (G.10), shows that the 
following relations exist between the physical and transformed spatial derivatives: 


(G.15) 

W (x 5 r i- r i;V 

(G.16) 


(G.17) 

’'x- -V ( V<r r 6V 

(G.18) 

’’r' VV'rW 


The denominator in Equations (G. 15) to (G.18) represents the area of the physical space in 
relation to the transformed space and is termed the inverse Jacobian, J -1 .where the direct 
Jacobian, J, represents the area of transformed space with respect to the physical coordi- 
nate system. Therefore, using the Jacobian, J, to represent (x^r^ - r^x^) -1 . Equation 
(G.7) becomes: 

(G.19) i(rQ) + r„ji(rE-rE v ) - r 5 jl(rE-rE v ) - x„J^(rF-rF v ) 

+ x 5 J^(rF-rF v ) + H - H v = 0 

Equation G.19 can be cast in conservation form as follows. First, note that the Jaco- 
bian is independent of time, therefore, (1/J) 3rQ/ (dx) = d (rQ/J) . Next, note that the 
spatial derivatives have the following characteristics: 
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(G.20) 


3 [rE - rE y ] 9[rE-rE y ] 
r TI 51 r l 3rj 


3 [r (rE-rE y )] 3 [r, (rE - rE y ) ] a 0 

31 5^ [rE_rE v ] ( 5! r iT3riV 


(G.21) 


3[rF-rF y ] 3[rF-rF y ] 

- x ri 51 + x i — ^ — 

3[x (rF-rF y )] 9 [x* (rF - rF y )] a a 

31 + 3rj + [TF - r ^ ] ^\-^ 


The last terms on the right-hand side of Equations (G.20) and (G.21) cancel out. There- 
fore, Equation (G.19) can be written as: 


(G.22) 



3[ rT1 (rE-rE v )] 3 [r^ (rE - rE y ) ] 


dx 31 3q 

3[x (rF-rF v )] 3 [x £ (rF - rF y ) ] H-H 

! 4- + = o 

91 3q J U 


In the radial momentum equation, the pressure gradient term can be written in a from sim- 
ilar to the fully conservative counterparts: 


(G.23) 


r r t 3 p 3P-] 

j[^ai + ^an. 



a(x T! P) 

“"ai - 


a(x^p)- 

arj 


Substituting Equations (G.15) to (G.18) into Equation (G.22) yields the transformed 
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Navier-Stokes equations: 


(G.24) 




K L 


^ x (rE-rE v )+^(rF-rF v ) 


(rE - rE v ) + T| r (rF - rF y ) 


dn 


H-H 

\ 

r~ 


= 0 


Equation (G.24) is the generalized form of the axisymmetric Navier-Stokes equations 
modelled in this study. Using over-script nomenclature, Equation (G.24) can be written in 
the following condensed form: 


(G.25) 


3(rQ) 

dx 


9[r(E-E v )] 5[r(F-F v )] 

51 + 5rj 


+ H - H v = 0 


where the overscripted terms are defined as follows: 


(G.26) 


(G.27) 


(G.28) 

t, 

E,= , 

(G.29) 

’lx E + 1r F 
F= J 

(G.30) 

n x E v + l1 r F v 
F v = , 

(G.31) 


(G.32) 

h v 3 
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The modelling of the viscous fluxes was performed in nonconservative terms, leaving 
the conservation aspect of the viscous terms to the proper evaluation of the metrics and 
flow difference procedures, as outlined in Section 3 and in Appendix B. For instance, the 
radial heat flux was represented by the following expression: 


(G.33) 
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A ppendix H 

Reference Plane Characteristic Method 

The method of characteristics has been successfully applied to numerous computa- 
tional techniques for modelling both the flow in the interior and at the boundaries. In this 
study, a reference plane characteristic method was applied at computational boundaries. A 
description of the reference plane method and its derivation from the governing fluid flow 
equations follows. 

In primitive form, the equations of mass, momentum, and energy (for axisymmetric 
flow) may be written as follows: 


(H.l) 


(H.2) 


(H.3) 


dp 

du 

dv 

dp 

dp pv 

= 0 

dt 

+p s 

+ p aF 

+ U^r— + V 

dx 

dr r 

du 

du 

du 

1 dP 

Viscous 

terms 

dt 

+ U 5J 

+ V dr 

+ 
-O 1 

ai 

ii 

P 


dv 

dv 

i n 

dv 

1 \j 

1 dP 

I — 

Viscous 

terms 

dt 

+ U d^ 

dr 

p dr 

P 



dP dP dP 2 

(H.4) -=r- +U^~ + V^r- - a 

dt dx dr 


'dp dp dp 
+ u^- +v-=fr 
dt dx dr 


= Viscous and source terms 


In generalized coordinates, x, and T[, Equations (H.l) to (H.4) become: 


(H.5) 8p 


Tz +<> 

du 


' du dv du dv' 

.^^ + ^ + ^ x 5q +T1 r ^_ 


+ u 


. 3p 3p 

S&, +11 x5n. 


+ V 


\ 9 P dp 

.^55 +11 rSn 


r 


ru*\ du r du du~l du du' 

(tl.O) — + u t ct + Tj =r- +V t =-=- +T1 =— 
dx 'xdriJ ’rdE, l rarj 


du du 


1 

+ - 
P 


' dP dP 

,^d^ +T1 x^j 


Viscous terms 


dv 

(H.7) s + u 


dv dv 

,^ +Tl x^ 


+ v 


dv dv' 

L^d? + r, rdn. 


i 

+ - 
p 


r dP dP 

3 r d? +Tlr dfl. 


Viscous terms 




(H.12) 

3P + U 3P - a 2 r 3p +U 3p l - V 
55 a 55J _ 4 

where the terms on the right-hand side of the above expressions are shown in Equations 
(H.13) to (H.16): 

(H.13) 

pv du d\ 9p 

l - — P^r^n V drj 

(H.14) 

w _ ^x3P Viscous terms 

2 p 8 t| dq + p 

(H.15) 

yjj _ ^r 5P ^ Viscous terms 

3 p dq dr| + p 

(H.16) ^ = 

n 

3P 9 3p ■v'' 

- V^ + a z V ^ + 2 ^ [y^T _ (Yj - 1 ) hj] 0 . + Viscous terms 


i = 1 


where V = uri x + vr] r is the contravariant velocity in the r|-direction. 
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A linear combination of Equations (H.9) to (H.12) can be formed by multiplying by 
the arbitrary weighting factors, /• (i = 1, 2, 3,4), respectively, and then summing. The 
resulting equation is: 


(H.17) / 1 


dp t dp du dv 

5f +u 35 +p ^g5 + P^5" ¥ i 


+ In 




+ L 




+ h 




= 0 


Equation (H.17) can be rearranged and put into a form that represents coefficients acting on 
the flow derivatives: 


(Zj - a 2 / 4 ) ^ + (U - a 2 U/ 4 ) | + Z 2 | + (P^ j + UZ 2 ) | 

+ 1R + Vi + U/ 3^ + l A + [ (^x /p ) *2 + (^r /p ^ h + U/ 4 ^ 




The following set of characteristic vectors is defined, which are the coefficients of the 
partial derivatives in Equation (H.18): 


(H.19) 

Wi = 

(U/, -a 2 U/ 4 -a 2 / 4 ) 

(H.20) 

w 2 = 

(py 1+ u/ 2 ,(j) 

(H.21) 

w 3 = 

(pVi+u/ 3 ,( 3 ) 

(H.22) 

W 4 = 

(^/ 2 /p+y 3 /p+u/ 4 ./ 4 ) 


Then Equation (H.18) can be written in the following form: 
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(H.23) dp| Wi +du| W2 + dv| W3 + dP| W4 = / 1 V 1 +/ 2 V 2 + / 3 V3 + / 4 ’P 4 

The characteristic normal, N, in the ^ - x plane is found by selection of the 
/j (i = 1, 2, 3, 4) such that the characteristic vector, Wj (j = 1, 2, 3, 4) , is normal to N. 
Mathematically, this condition requires that the dot product of the characteristic normal 
and characteristic vector be equal to zero. 


(H.24) 


N • Wj = 0 (j = 1, 2, 3, 4) 


In matrix form. Equation (H.24) is: 


(H.25) 


UN$+N t 

0 

o -a 2 (UN^+N T ) 

V 


UN^+N t 

0 0 

h 

p5|N$ 

0 

UN^+N t 0 

h 

0 

Wp 

Wp UN^+N t 

J 4_ 


For a nontrivial solution, the determinant of the coefficient matrix must be equal to 
zero, leading to the following requirement: 


2 

(H.26) ( UNj= + N t ) [(UN^ + N x ) 2 -a 2 ^ N 2_ a 2^2 N 2 j = Q 

Equation (H.26) requires that one of the following two expressions must be true: 

(H.27) UNj: + N t = 0 

(H.28) UN^ + N t = ± aN^ + £2 

Noting that d^/dx = - N x /N^ , Equations (H.27) and (H.28) can be written as: 



(H.29) 


d| 

dx 


d| 

dx 


U 


(H.30) 


U ± 


Equation (H.29) represents the projection of the pathline onto the r| = constant planes, and 
Equation (H.30) represents the projection of the Mach cone onto the T| = constant plane. 
Substitution of Equation (H.29) into Equation (H.25) yields the following expression: 


(H.31) 


- 



- 



0 

0 

0 

0 


/, 

P^x N 5 

0 

0 

0 


h 

P^r N q 

0 

0 

0 


h 

0 

Wp 

Wp 

0 


U 


Solving Equation (H.31) for l x to / 4 gives: 

(H.32) 1 1=0, — - — / 2 h — — l ^ = 0, / 4 is arbitrary 

Two possible solutions are l x = / 4 = 0 and l 2 arbitrary, or /, = / 2 = / 3 = 0 and / 4 arbitrary. Sub- 
stituting /j = / 4 = 0, arbitrary, and / 3 = - ^ x / 2 /^ r , into Equation (H. 17) yields one of 
the two compatibility equations that are valid along the characteristic curve 
dtjdx = U: 


(H.33) 


t £ t £ 

. du ... du ^x dv ... ^xdv 1 dP 1 dP 

2 dx + 2 d£ 2 dx 2 £^d£, + p 2 d£, p 2 d^ 


^x 

^r 
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Dividing through by the 
(H.34) gives: 

nonzero arbitrary multiplier l 2 and then simplifying Equation 

(H.34) 

£ du-£ dv = (^ 'P, 'P,)dx 

where 

• 

(H.35) 

, du du di: 

du = y dx + vr — — dx 
dx dq dx 

(H.36) 

dv dv d£ 

dv — y dx + »T7 -t— dx 

dx dq dx 


Substituting /] = / 2 = / 3 = 0 and l 4 arbitrary into Equation (H.17) yields the other com- 
patibility equation that is valid along the characteristic curve d£/dx = U. 


(H.37) 


_ _ a 2 U^B + — + U— - 

a ax a u aq + ax ^ ^ 




'4 = 0 


Dividing through by the nonzero arbitrary multiplier / 4 and putting the source term on the 
right-hand side. Equation (H.37) gives the compatibility equation: 


(H.38) 


^dp 9 3p dP dP 

“ a ai “ a u a? + ai + u a? = ^4 


Along the characteristic line d^/dx = U , Equation (H.38) can be expressed in terms of 
total derivatives: 


(H.39) 


dP - a 2 dp = ^ 4 dx 


Two additional compatibility equations are obtained for the characteristic curves spec- 
ified by Equation (H.30). Substitution of Equation (H.30) into Equation (H.25) yields: 



0 


0 


(H.40) 


P^x N 5 ±aN^ + ^ 

o 

0 ^x N 5 7 P 


±aNt 



0 


Wp ±aN Ji 


Solving Equation (H.40) yields: 


(H.41) 



/ = - p ^' = pa ^' j 

3 ' 

l 4 is arbitrary 


Substituting Equation (H.41) into Equation (H.17) yields: 
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(H.42) 


2t , P aU ^x^4^9u 

'JR* (’ 


P a V4 3v / 2 f , P a ^V4\dv 


P aU V4\dv 


. ap 

+ , 4H + 


l-Uf 


<' 4 1 


£F«?. 


ap 


pa 5/4 P a V 4 

a 2 / 'P. T - — — V, T - 7 

4 ' ^ 


'VV^ 


Substitution for the total derivatives, such as dP/dx = (dP/31;) d^/dx + 9P/dx , 
where d^/dx = U qF a J^ 2 + £ 2 , and likewise for du and dv. Equation (H.42) simplifies 
to: 


(H.43) dP =f 




a 2x Pj t 



T 



u/ ,u> 

3 4 


which is the compatibility equation that applies along the projection of the Mach cone on 
the constant-T] reference plane. 

Summarizing, the applicable characteristic and corresponding compatibility equa- 
tions valid on a constant-T) reference plane are given in Table 26. 
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Table 26. Characteristic and compatibility equations for a constant-q reference plane. 


Characteristic 

equation 

Compatibility equation 

*-u 

dx 

^ 

£-u 

dx 

dP - a 2 dp = *P 4 dx 

§ “ V ~ a fcl + % 

pa£ x pa^ r 

dP - , 1_ du dv 

= fa 2 '*' - pa ^ x v 2 - ^3 + sOdx 

l &+% £ + s? v 

1 « + 

P a £ x pa^ r 

dP + dU + :"-v - dV 

- (a 2 'i‘. + Pa S .¥, + f ^_'P 3 + 'r 4 ldt 

l ’ £♦«? 2 £ + «? 3 4 J 
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The right-hand side terms, to T 4 ? are listed for completeness. 


(H.44) 

pv du dv 8p 

1 r aq 

(H.45) 

^ _ A ,au ^l x ap Viscous terms 
2 aq p aq + p 

(H.46) 

w - \r^ dP Viscous terms 

3 an p an + p 

(H.47) 

ap 9 ap 

*F. = - V^r- + a 2 V^ 

4 an an 


n 

+ [YjRjT - (y i — 1 ) hj] Oj + Viscous terms 

i = 1 

In all cases studied, the right-hand side terms were numerically approximated by a 
central difference approximation and were averaged between the boundary point and the 
point adjacent to the boundary. Metric differences were always evaluated within the com- 
putational region used by the dependent variables. Since the terms in Equations (H.44) to 
(H.47) were central differenced, the option of applying artificial diffusion was employed 
by adding the following expressions to the right-hand sides of the respective equations: 

(h.48) 

where T is given by: 


IVI+a^ . |U| + a 
Aq + AE, 


(H.49) 
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with |U| and | V| being the absolute values of contravariant velocity in the £ and q direc- 
tions, respectively. The parameter ns an artificial viscosity coefficient equal to the 
value applied in the free stream. The specific terms upon which the artificial diffusion 
expressions act (given by the Q symbol in Equation (H.48)) are listed in Table 27. 


Table 27. Artificial diffusion applied to constant-ri reference plane characteristic 

equations. 


Right-hand side term 

Term upon which 

♦ 

artificial diffusion acts 


P 

1 

^2 

u 

1 

¥ 3 

V 

1 

¥4 

p 

1 

p 

- a 2 


The viscous terms in the constant-q reference plane characteristic method simplify to 
the following terms for the axial momentum, radial momentum, and the energy equations. 

... „ q 9x q 9x x £ dx £ dx 

Axial momentum _ Jr xr + Jx xx ( xr ; Jr xr Jx xx 

viscous terms P P ^ pr P ^ p d£ 

n q 3x , q dx x £ dx £ dx x flfl 

Radial momentum _ Jx xr Jr rr + _rr Jx xr + Jr rr _ _J0 

viscous terms P P dq pr p d£ p d£ pr 


(H.50) 

(H.51) 
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(H.52) 


Energy transfer viscous terms 


dv 

8u 

f 8u 

dv V 

X T| •= r — 

it 'r an 
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c 8u 
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xr 
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^ , 3 <*x , 



The corresponding equations for the shear stress and heat flux terms are presented in 
Section 2. The evaluation of terms containing only T| derivatives was performed by aver- 
aging at the points on the boundary and the £ point closest to the boundary. The evalua- 
tion of the mixed derivatives was performed at the point midway between the boundary 
point and the next nearest interior £, point, whereas the evaluation of terms containing 
only E; derivatives was performed about the next nearest interior £, point. No effort was 
made to make the viscous terms one-dimensionally conservative. However, the fundamen- 
tal rule of viscous flow conservation, that metric grid points outside of the computational 
domain not be used, was adhered to in all cases except those terms containing second- 
order derivatives in the ^ direction. 

When applied along a solid wall or a centerline, the viscous terms and the species 
source function terms were assumed to be negligible. 

The characteristic and compatibility equations applicable to the steady state, inviscid, 
species equations are derived in Reference 23. Simple inspection allowed the extension of 
those results to transient, diffusive flow. The resultant characteristic and compatibility 
equations are given in Table 28. 
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Table 28. Characteristic and compatibility equations for species equations. 


Characteristic equation 

Compatibility equation 

« = u 

dz 

dCD i = V 4 + i dx 


The right-hand side term, ¥ 4 + j , is given below. 

9co. 

e 

35 


(H.53) ¥ 4 + i " £[ " P v ^ - ~ ^ 


*^ij jij 

x SFff " n r 35j " T 


The equations for the diffusive mass flux terms are listed in Section 2. 

The handling of the advective and diffusive terms was performed in a manner analo- 
gous to the terms arising in the mass, momentum, and energy equations. Likewise, an arti- 
ficial diffusion option exists for the species transport equations, being identical in form to 
the equations for mass, momentum, and energy. Table 29 lists the variables to the species 
equations upon which the artificial diffusion acts. 

Table 29. Artificial diffusion applied to the species equations. 


Right-hand 
side term 

Term upon which 
artificial diffusion acts 


V 4 + i 

(D. 

l 

1 


The analogous constant- reference plane formulation is readily derived from the 
constant- rj reference plane formulation derived in this appendix, by simply interchanging 
the \ and T| coordinates wherever they appear in all of the equations. 












237 


A ppendix I 

Chemical Reaction Mechanism and Reaction Rate Constants 

Eight possible chemical reactions were considered as potentially significant to the 
hydrogen-oxygen chemical rocket engine. All of the backward chemical reaction rates 
were modelled as a modified Arrhenius equation: 

(1.1) K bj = BjT^exp 

where j denotes the specific chemical reaction, and Bj, (Xj, and Ej are determined from 
experimental data. The plausible chemical reactions and reaction rate constants employed 
in the Tethys code (i.e., the default values) are listed in Table 30. The values of the reac- 
tion rate constants were obtained from recommendations made by Bittker (Reference 55) 
in 1976. The units in the following table are presented in terms of cubic centimeters, 
degrees Kelvin, calories, moles, and seconds. 


Table 30. Chemical reaction rate constants. 


Chemical Reaction 

B J 

a. 

j 

E i 

h + h + m<=*h 2 + m 

2.2 x 10 14 

0.0 

96,000.0 

H + OH + M <=> H 2 0 + M 

1.3 x 10 15 

0.0 

105,140.0 

0 + 0 + M <=> 0 2 + M 

1.8 x 10 18 

- 1.0 

118,000.0 

OH + M<=»0 + H + M 

7.1xl0 18 

- 1.0 

0.0 

o+oh<=*o 2 + h 

1.89 x 10 14 

0.0 

16,400.0 

H + OH <=>H 2 + 0 

4.2xl0 14 

0.0 

13,750.0 

H 2 0 + HoH 2 + 0H 

4.74X10 13 

0.0 

6098.0 

oh + oh»h 2 o+o 

6.8xl0 13 

0.0 

18,365.0 
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The default values of the catalytic efficiencies (i.e., B i j / B j , as discussed in Section 
2.3) used in the Tethys code were obtained from the same reference as the chemical reac- 
tion rate constants (Bittker, Reference 55). They are shown in Table 31. 

Table 31. Catalytic efficiencies. 


Chemical Reaction 

H 

h 2 

h 2 o 

O 

o 2 

OH 

h + h+m«=>h 2 + m 

1.0 

4.1 

15.0 

1.0 

2.0 

1.0 

h + oh+m<=>h 2 o + m 

1.0 

4.0 

20.0 

1.0 

1.5 

1.0 

O + O + M <=>0 2 + M 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 

OH + M<=>H + 0 + M 

1.0 

1.0 

1.0 

1.0 

1.0 

1.0 
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A ppendix J 

Chemical Kinetics Validation 


To evaluate the accuracy and methodology employed in modelling the chemical spe- 
cies source term, a comparison was made between a well-established one-dimensional 
chemical kinetics model, LSENS (Reference 19), and the model employed in Tethys. The 
geometry and boundary conditions considered in the flow verification study correspond to 
a supersonic source flow with an area ratio of 100 and a total length of 1 foot, as shown in 
Figure 84. The inlet boundary conditions employed for the verification roughly corre- 
spond to the equilibrium conditions that would exist near the throat of a hydrogen-oxygen 
rocket engine operating at a H 2 /02 mixture ratio of 8.0 and a fluid static density near 0.05 
lbm/ft 3 . The initial velocity was specified such that the initial conditions corresponded to 
an inlet Mach number of 1.5. Three different inlet fluid densities, encompassing two 
orders of magnitude in difference, were studied to assure that a rather complete spectrum 
of chemical recombination rates and mechanisms could be studied. The inlet conditions 
are presented in Table 32. Reference plane characteristics, with chemical kinetics, were 
used at the exit plane. 



r i 



Figure 84. Schematic of source flow geometry used for chemical kinetics validation. 
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Table 32. Inlet boundary conditions for reacting source flow study. 


Property 

Inlet value 

Units 

T 

6060.6 

R 

P 

0.005, 0.05, 0.5 

Ibm/ft 3 

V 

7171.9 

ft/sec 

% 2 

0.0156 

Dimensionless 

“V 

0.7749 

Dimensionless 

“o 2 

0.0829 

Dimensionless 

®h 

0.0025 

Dimensionless 

“oh 

0.1050 

Dimensionless 

®0 

0.0191 

Dimensionless 


In all of the cases examined, a 21 X 4 grid point structure was employed in the 
Tethys analysis, as shown in Figure 85. In general, the one-dimensional flow model, 
LSENS, used about 500 marching steps. Due to the much larger number of steps taken and 
the maturity and higher-order accuracy of the LSENS code, it will be assumed to be 
numerically exact, with the understanding that the coefficients used for the actual chemi- 
cal reaction rates could be significantly different from those that exist in reality. 




Height, ft 
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Radius, ft 


Figure 85. Computational grid used for the chemical kinetics source flow validation. 


c;r> 
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In Figure 86, the variation of Mach number with the radius ratio (or equivalently, the 
area ratio), as predicted by the LSENS flow model, is shown for the three different inlet 
flow cases examined. It is readily observed that the differences in the Mach numbers are 
fairly large, and solely the result of differences in the chemical reaction mechanisms 
occurring at the different fluid densities. The tendency for increased chemical recombina- 
tion to occur at higher densities results in an effective energy release, very similar to a 
Rayleigh Line effect, with a resultant lowering in fluid Mach number with increasing 
recombination. 



Low density 

Moderate density 
High density 


Figure 86. Mach number distribution from LSENS code for chemical kinetics study. 


Figure 87 presents the Mach number error associated with using the 21 point grid in 
the Tethys code. The result are quite good, with the maximum error at any location being 
less than 3 %. The error is less than 1 % at the nozzle exit. 
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The results for the high density case suggest that somewhat greater packing of the 
grid points near the inlet would be beneficial. In fact, approximately twenty percent of the 
total number of steps taken by the LSENS code were taken between the inlet and the first 
point used in the Tethys code. It is evident that the species chemical kinetics reaction rates 
and source term evaluation applied in the Tethys code is very good at predicting the flow 
Mach number and would correspondingly be expected to accurately predict other fluid 
properties of interest such as temperature and pressure. 

In Tables 33 to 35, the chemical species generation rate at the nozzle inlet is presented 
for the Tethys and LSENS codes. It is evident that the species generation rates are essen- 
tially the same, with the slight differences in generation rates (< 2 %) possibly being due 
to round-off errors associated with converting the chemical reaction rates and chemical 
equilibrium constants to consistent units. 
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Table 33. Species generation rates at the inlet, low density analysis. 


Species 

Cj, LSENS code, lbm/(ft 3 -sec) 

Oj, Tethys code, lbm/Cft^-sec) 

h 2 

-2.6694X10 1 

-2.6666X10 1 

°2 

3.2694X10 1 

3.2643X10 1 

OH 

1.8899X10 1 

1.8369X10 1 

H 

1.7509X10 1 

1.7463 X10 1 

0 

- 1.1449X10 2 

- 1.1433X10 2 

h 2 o 

7.2087 X10 1 

7.2521 X10 1 


Table 34. Species generation rates at the inlet, moderate density analysis. 


Species 

a,, LSENS code, lbm/(ft 3 -sec) 

Gj, Tethys code, lbm/(ft 3 -sec) 

h 2 

-2.6214X10 3 

-2.6190X10 3 

°2 

3.2716X10 3 

3.2665 X10 3 

OH 

- 2.4478 X10 3 

- 2.4588 X10 3 

H 

1.4230 X10 3 

1.4213X10 3 

O 

- 1.1633 X 10 4 

-1.1617 X 10 4 

h 2 o 

1.2008 X10 4 

1.2007 X10 4 
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Table 35. Species generation rates at the inlet, high density analysis. 


Species 

Oj, LSENS code, lbm/(ft 3 -sec) 

Oj, Tethys code, lbm/(ft 3 -sec) 

h 2 

-2.1414X10 5 

—2.1431 X 10 5 

o 2 

3.2939X10 5 

3.2885 X10 5 

OH 

-4.5825 X10 6 

-4.5416X10 6 

H 

-1.8562X10 5 

-1.8290X10 5 

0 

-1.3468X10 6 

- 1.3451 X10 6 

h 2 o 

5.9996X10 6 

5.955 1X10 6 


In Figures 88 to 90, the radial variation (as predicted by the LSENS code) of the radi- 
cals OH, H, and O, are shown as the flow expands to an area ratio of 100, for the three dif- 
ferent inlet fluid densities studied. From these figures it is readily apparent that the fluid 
density has a profound effect on the concentration of the radicals, especially monatomic 
hydrogen and monatomic oxygen. 
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In Figures 91 to 93, the error associated with using the Tethys code with the 21 point 
computational grid is presented. In these figures it is seen that the errors in predicted radi- 
cal species concentrations is of the same general order of magnitude as the errors in pre- 
dicting the Mach number presented in Figure 87, with the maximum errors in species mass 
fraction being less or about 5 % and the maximum error for any given species always 
occurring with the moderate density test case. 

In view of the rather accurate ability of the Tethys model to predict both the funda- 
mental flow properties and radical species concentrations, even with a rather crude com- 
putational grid, it is concluded that the species source term is well modelled in the Tethys 
code. 



Figure 91. Error in OH mass fraction in a chemically reacting source flow study. 
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Figure 92. Error in monatomic hydrogen mass fraction in a chemically reacting 
source flow study. 



Figure 93. Error in monatomic oxygen species mass fraction in a chemically 
reacting source flow study. 
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The verification of the fluid properties methodology used in the Tethys model is 
essential to the accurate prediction of rocket engine performance and heat transfer. Some 
of the properties, such as specific heat, viscosity, and thermal conductivity, and their corre- 
sponding mixture expressions, can be readily verified with existing and mature computer 
models which employ the same expressions. On the other hand, other coefficients such as 
the multispecies diffusion coefficient, the thermal diffusion coefficient, and the Dufour 
coefficient, are much more difficult to verify, requiring one to rely on sparse data or engi- 
neering judgement. 


K.l. Validation of Specific Heat, Entropy, Thermal Conductivity, and Viscosity 
The equations governing the definition of the multispecies specific heat and entropy 
are based on thermodynamic arguments and are discussed in Section 2 and Appendix B, 
respectively. The multispecies mixture equation employed for the viscosity and thermal 
conductivity is of an approximate nature and is given below: 


(K.l. 1) 



where the denominator, |3., for the fluid viscosity is defined as: 


(K.1.2) Pi=V £ VP* 

j= 1 
j*i 

and the empirical expression for the viscosity interaction coefficient, (p^, is given below: 
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The mixture formula used for the thermal conductivity is given in Equation (K.1.4): 


n 


(K.1.4) 


Xjkj 

“St 


where the denominator, fy, for the mixture thermal conductivity is defined as: 


(K.1.5) Pi=V 

j= 1 
j*i 

and the empirical expression for the thermal conductivity interaction coefficient, 'E.j, is 
given by Equation (K.1.6): 


(K.1.6) 



where (p- is the viscosity interaction coefficient described in Equation (K.1.3). 

To validate that the Tethys model employs proper techniques for predicting the fluid 
specific heat, entropy, thermal conductivity, and viscosity, two computational cases were 
studied. In one case, the fluid was near the combustion temperature of a hydrogen-oxygen 
rocket engine (7000 R). In a second case, the fluid was at 1710.5 R. In both cases, the fluid 
had an oxygen to fuel ratio of 8.0 and the species concentration was in chemical equilib- 
rium. These two cases represent conditions that are at the high-temperature and the low- 
temperature portions of the thermodynamic curve fit, with the split occurring at 1800 R. 
The inlet fluid conditions are given in Table 36. Table 37 compares the fluid properties 
prediction from CET89 (Reference 47) and the Tethys code. The values are seen to be in 
excellent agreement, with any slight discrepancies undoubtedly due to the fact that the 
thermodynamic data file used in the Tethys code was an older version of the data file used 
by CET89. 
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Table 36. Fluid conditions for validation of specific heat, thermal conductivity, and 

viscosity models. 


Fluid property 

High temperature 

Low temperature 

Temperature, R 

7000.0 

1710.5 

X 

3 

0.01507 

0.02600 

8 

o 

to 

0.07910 

0.12281 

<°h 2 o 

0.78276 

0.58549 

8 

O 

SC 

0.10308 

0.19312 

C0 H 

0.00235 

0.00815 

<°0 

0.01764 

0.06443 


Table 37. Comparison of CET89 and Tethys fluid property predictions. 


Fluid property predicted 

High temperature 

Low temperature 

CET89 

Tethys 

CET89 

Tethys 

Specific heat, Cp,Btu/lbm-R 

0.7819 

0.7702 

0.5570 

0.5572 

Entropy, s, Btu/lbm- R 

4.285 

4.274 

3.865 

3.863 

Thermal conductivity, k, Btu/ft-sec X10 5 

9.367 

9.349 

1.841 

1.839 

Viscosity, |i, lbm/ft-sec X10 5 

7.521 

7.520 

2.448 

2.448 


K.2. Validation of Multiple Species Diffusion Coefficients 
Lacking an established computer model upon which to verify the accuracy of the 
multiple species diffusion coefficient model, the Tethys model was verified by comparing 
the computational results of a ternary system with the theoretical results. Since the ternary 
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case represents a multispecies condition that computationally differs from a six species 
case only in the number of computational loops required in the Tethys code, validation of 
this case is a reasonable validation of the entire multispecies calculation procedure 
involved in Tethys. 

The exact formulation of the multispecies diffusion coefficients for a ternary system 
is extremely lengthy, being presented in Reference 9. The ternary system considered com- 
prised hydrogen, oxygen, and steam at a temperature of 1000 R. The values of the binary 
diffusion coefficients used in the theoretical results and the Tethys model were obtained 
from curves fits of tabulated values suggested by Svehla (Reference 56), who employed a 
combination of theory and experimental data. The results of the study are presented in the 
Table 38 and are shown to be in excellent agreement. 


Table 38. Comparison of theoretical and predicted multispecies diffusion coefficients for 

a ternary system. 


Multispecies diffusion 
coefficient 

Theoretical value 

Tethys prediction 

A H 2 -H 2 

-1.443 x 1(T 5 

-1.443 x 10 -5 

a h 2 - 0 2 

-5.350 x 10 -5 

-5.351 x 10 -5 


K.3. Validation of Thermal Diffusion Coefficient Procedure 
Verification of the accuracy of the thermal diffusion model used in the Tethys code, 
along with the accuracy to which the multiple equations were solved was demonstrated in 
a rather heuristic manner. First, the case for which binary thermal diffusion data exists, 
that of the H 2 /O 2 species system, was examined and compared with the predictions from 
the Tethys model. Second, an H 2 /O 2 system, complete with reacted and dissociated spe- 
cies present, was examined at numerous temperatures for any irregularities in the shape or 
magnitude of the diffusion coefficient curves. 
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In Figure 94, the predicted variation of the thermal diffusion coefficient for an oxy- 
gen/hydrogen system from a temperature range of 160 to 530 R is presented, which corre- 
sponds to the temperature at which the mean thermal diffusion factor was .measured. The 
integrated average value of the thermal diffusion coefficient predicted by Tethys is 0.206 
while the measured value reported in Reference 57 is 0.192. The agreement between pre- 
dictions and experimental measures is therefore taken to be very good in view of the com- 
plexity of the driving forces describing the nature of the thermal diffusion coefficient. 



Figure 94. Comparison of experimental and predicted values of the thermal 
diffusion factor in a hydrogen/oxygen system. 

A final heuristic validation of the thermal diffusion modelling is presented in Figure 
95, where an oxygen/hydrogen system (at the stoichiometric mixture ratio of 8.0), of the 
same chemical composition as tested in the chemical kinetics validation studies (Appen- 
dix J), is presented. Several conclusions can be drawn from this figure. The first is that the 
curves appear to be regular in nature and are about an order of magnitude less than that of 
an unreacted hydrogen/oxygen system. This is as expected since the mass fractions of 
most of the constituents in the reacted stoichiometric system are rather small. A final 
observation to be made is that the thermal diffusion coefficient curve for steam is negative 
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at lower temperatures and becomes positive at temperatures near 1000 R. This is expected 
since the binary diffusion theory would suggest that the water vapor diffusion coefficient 
would be a minimum near 700 R, which is in agreement with the predictions shown in 
Figure 95. 



O 2000 4000 6000 8000 

Temperature, R 

Figure 95. Temperature sensitivity of the thermal diffusion coefficient 
in a chemically reacted hydrogen/oxygen system. 
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To validate the accuracy of the Tethys model in describing the transfer of mass as the 
result of ordinary diffusion, a simple flow model was considered. The viscous flow of a 
fluid in a pipe, subjected to slip boundary conditions along the wall and a trace concentra- 
tion of a species at the wall was modelled. This case is entirely analogous to the case of a 
solid bar at a prescribed temperature that is suddenly subjected to a different temperature 
at the outer surface. The corresponding computational grid and boundary conditions 
imposed for this test case are presented in Figure 96 and Table 39. 


o 

§1 



o. o.i 

Axial distance, ft 


Figure 96. Computational grid and geometry of species diffusion validation study. 



256 


Table 39. Boundary conditions for species ordinary diffusion validation study. 


Inlet 

Exit 

Centerline 

Wall 

p = 0.01 lbm/ft^ 

dP - o 
35~ 0 

0 

II 

Oh I P" 
ro |rt> 

3p _ 0 

5r| “ 0 

u = 1000 ft/ sec 

O 

II 

0 

II 

|ro 

II 

O 

v = 0 

O 

II 

v = 0 

v = 0 

T = 500 R 

0 

II 

fel# 

0 

II 

USs” 

II 

0 

(D„ = l.Oe - 06 
”2 

d \, 

K 2 “° 

3 “h, 

^r- =0 

(Ox, = l.Oe — 04 
H 2 


The results of the Tethys code and theoretical results (extracted from graphical data. 
Reference 13) are shown in Figure 97. The results at distances greater than 0.025 ft are in 
excellent agreement with theory, with the most likely reason for the discrepancy at x = 
0.025 ft being that the graph from which the theoretical predictions were made was based 
on the assumption that the Fourier number was rather large. 
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The prediction of rocket nozzle performance parameters was performed in a manner 
which conformed with the JANNAF methodology (Reference 58) for nozzle performance 
prediction as closely as possible. The procedure consists of using a control surface that 
crosses the throat of the nozzle and envelopes the divergent section of the nozzle and the 
complete combustor and inlet section, as illustrated in Figure 98. 



Figure 98. Control surface employed for evaluation of nozzle performance parameters. 

Using the control surface shown in Figure 98, the prediction of nozzle performance is 
obtained by evaluating the pressure forces, momentum flux, and mass flux that cross the 
throat and the pressure forces and shear forces that act on the divergent section of the noz- 
zle. The procedure for the evaluation of these terms follows. 

M.l. Evaluation of Nozzle Throat Region Parameters 
Consider a control surface cutting across the throat region of a rocket nozzle, as illus- 
trated in Figure 99, where a constant-^ line, originating at point o on the centerline, strikes 
the nozzle wall at the throat, point e. 
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I 

I 

o 




Figure 99. Control surface across the nozzle throat. 


The mass flow rate across control surface oe is given by the following equation: 


(M.1.1) m = Jp (V • n) dA = Jpu7td(r 2 ) - J pv27trdx 


To maintain consistency with the manner in which the finite difference equations were 
applied, the numerical evaluation of the mass flux was performed by a centered integration 
technique, as shown in Equation (M.1.2), where the total number of computational points 
located on a constant £ line is denoted by N2, with individual points on that line described 
by the subscript j. 


(M.1.2) 


N2- 1 


m«JC ^jT (pu) | 
j = 2 


r lj+l- r lj-l 


N2- 1 
-2lt £ 

j = 2 


(pvr) |. 


xl j + l 


- x 


j-1 


JL 


Edge 

points 


The axial momentum flux crossing the control surface can be estimated in a manner 
analogous to that in which the mass flow rate is calculated. The resulting expression is 
given by Equation (M.1.3): 



260 


(M.1.3) 


Axial 

momentum ~ k 
flux 


N2- 1 


X 


r_2| 


[ (pu)u] 


'j+l 


- r 


‘j-l 


N2- 1 

_ 27t X t (p vr > ij 

j = 2 


xl j+ 1 ~ xlj_ i 


Edge 

points 


In a similar manner, the contribution of the axial pressure force component acting on 
the throat control surface can be approximated by Equation (M. 1.4): 


(M.1.4) 


Axial 

pressure 

force 


r e N2- 1 

= JPttdr 2 = n Plj 

r o j = 2 




+ 


Edge 

points 


In Equations (M.1.2) to (M.1.4), the edge points are all evaluated using a first-order 
method. 


M.2. Evaluation of Nozzle Divergent Wall Force Components 
In order to maintain consistency with the JANNAF nozzle evaluation procedure, it is 
necessary to evaluate the forces that act upon the nozzle divergent section. The forces are 
comprised of two component, a thrust contributing component from the fluid static pres- 
sure and a degradation component arising from the fluid shear stresses. 

The axial component of thrust obtained by the fluid static pressure forces acting on 
the nozzle surface can be derived in a manner that is entirely analogous to the procedure 
employed in establishing the pressure forces that act on the control surface crossing the 
throat. The resulting expression is given by Equation (M.2.1): 


*N1 


Ni-i r r 2 


= J P7trdr 2 = 7t Plj 


*T1 


i = T1 + 1 


i+1 r 1 i — 1 


Edge 

points 


(M.2.1) 


Axial 

pressure 

force 
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where N1 denotes the total number of computational points on the surface of the nozzle 
wall and T1 corresponds to the computational point that lies at the throat of the nozzle. 

The total amount of drag that acts in the divergent region of the nozzle is given by 
Equation (M.2.2): 


Axial 
(M.2.2) drag 
forces 


x = N1 

f x2ttrdx 
x =T1 


Nl- 1 


Y [x27tr] J 
i = T1 + 1 



+ 


Edge 

points 


The evaluation of the drag forces that act on the surface of the nozzle is somewhat 
more complicated than the expressions given for the pressure force terms since the shear 
stresses are actually calculated in the Tethys model at the point (see Appendix B) midway 
between the nozzle surface and the adjoining interior point, as illustrated in Figure 100. 




Nozzle surface 

Location at which 
shear stress is determined 


i — 1, j — 1 i, j - 1 i+l.j-l 


Figure 100. Nozzle surface examination of shear stress calculation procedure. 


Consequently, the procedure for determining the shear stresses that act on the nozzle wall 
is to examine the manner in which the averaged shear stress is calculated and to determine 
the relation of this term to that which exists along the surface of the nozzle wall. The 
dependence of the wall shear stress on a nearby point can be determined through consider- 
ation of a locally transformed grid. Assuming that the flow is parallel to the nozzle surface 
and that the difference in radius between the wall point and the nearest interior point is 
negligible, the transfer of momentum from the fluid to the wall can be determined by solu- 
tion of the following equation (Reference 59, page 474): 


(M.2.3) 



dP 

at 


where the t subscript denotes the direction tangent to the wall and the n subscript denotes 
the direction normal to the wall. Assuming that the tangential pressure gradient at the wall 



262 


is constant from the wall surface to the nearest midway interior point (an approximation 
consistent with the assumption of parallel flow). Equation (M.2.2) can be integrated and 
rearranged to yield the following: 


(M.2.4) 




iJ-1/2 


dP 

~5t 


An 


The evaluation of the pressure gradient term in Equation (M.2.4) is a relatively 
straightforward procedure, with the normal distance corresponding to 1/2 the distance 
from the wall point at the i, j position to the interior point located at the i, j - 1 location. 

The evaluation of the shear stress term located at the i, j - 1/2 point is performed in 
a manner consistent with the procedure used in Appendix B. Using this procedure, the 
shear stress terms at the midway point is equal to the following: 


(M.2.5) 


tn 



± ln 2 + n 2 
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i. j - 1 
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where 



(uri r + vTi x ) 



, is the velocity in the tangential direction. 



263 


The evaluation of the metric terms contained in Equation (M.2.5) obeys the funda- 
mental philosophy described in Appendix B which, in essence, means that all of the grid 
points used in the determination of the grid metrics arc never outside of the region in 
which the dependent variables are evaluated. 

In cases where Dirichlet boundary conditions were used for the thermal energy equa- 
tion, a procedure analogous to that used for the evaluation of the shear stress terms was 
employed for the evaluation of the heat transfer rates. 

M.3. Comparison with Experimental Results 

In the same manner in which all reasonable attempts to conform with the JANNAF 
nozzle evaluation procedure were applied to the Tethys code predictions, conformity of 
the experimental results with the prescribed JANNAF rocket nozzle test evaluation meth- 
odology (Reference 60) was also pursued. For the most part, the experimental results that 
were published were in agreement with the prescribed mannerism of publishing nozzle 
performance, although the need to maintain consistency between analytical and experi- 
mental results required some changes to the experimental performance measurements. 

The performance measurement alteration for the 1030:1 high area ratio nozzle tests 
was that the nozzle experimental evaluation program used an effective total pressure cor- 
rection concept that was inconsistent with the JANNAF methodology and was not applied 
in this study. 

The 25 lbf film cooled nozzle study, although not stated in the report, used chamber 
pressure measurements from a location just upstream of the convergent section to report 
characteristic velocity performance values. In this study, the chamber pressure measure- 
ments located at the inlet plane of the combustor were used to compute the nozzle stagna- 
tion pressure. 

No alteration of the heat flux measurements was made. In the high area ratio nozzle 
tests, the measurements were included in the report. In the 25 lbf film cooled nozzle study, 
the measurements were extracted by the author from unpublished coolant water mass flow 
rate measurements and the measured inlet and outlet temperatures of the water. 
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The methodology for applying artificial diffusion in the Tethys code is identical to 
that employed in the Proteus code (Reference 5), which is presented in this appendix and 
is a reproduction of Section 8 in Reference 5. The references cited within Section 8 of the 
Proteus manual are cited at the end of this appendix. Although nomenclature represent 
those employed in the Proteus report, they are consistent with those used throughout the 
present work, with only one exception. In the Proteus report, the nomenclature used for 
the spectral radii term is y, while in Appendix H the spectral radii is presented using the 
Greek character T. 
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8.0 ARTIFICIAL VISCOSITY 

With the numerical algorithm of Section 7.0, high frequency nonlinear instabilities can appear as the 
solution develops. For example, in high Reynolds number flows oscillations can result from the odd-even 
decoupling inherent in the use of second-order central differencing for the inviscid terms. In addition, 
physical phenomena such as shock waves can cause instabilities when they are captured by the finite dif- 
ference algorithm. Artificial viscosity, or smoothing, is normally added to the solution algorithm to suppress 
these high frequency instabilities. Two artificial viscosity models are currently available in the Proteus 
computer code - a constant coefficient model used by Steger (1978), and the nonlinear coefficient model 
of Jameson, Schmidt, and Turkel (1981). The implementation of these models in generalized 
nonorthogonal coordinates is described by Pulliam (1986b). 

8.1 CONSTANT COEFFICIENT ARTIFICIAL VISCOSITY 

The constant coefficient model uses a combination of explicit and implicit artificial viscosity. The 
standard explicit smoothing uses fourth-order differences, and damps the high frequency nonlinear insta- 
bilities. Second-order explicit smoothing, while not used by Steger or Pulliam, is also available in Proteus. 
It provides more smoothing than the fourth-order smoothing but introduces a larger error, and is therefore 
not used as often. The implicit smoothing is second order and is intended to extend the linear stability 
bound of the fourth-order explicit smoothing. 

The explicit artificial viscosity is implemented in the numerical algorithm by adding the following terms 
to the right hand side of equation (7.5a) (i.e., the source term for the first ADI sweep.) 

P (2) A _ ,W At . 

-Zj- (V{A { Q + V 7 A„Q) - -Zj- [(V^) 2 Q + CV/Q] (8.1) 

where e£> and ef 1 are the second- and fourth-order explicit artificial viscosity coefficients. The symbols V 
and A are backward and forward first difference operators. Thus, 

V {Q, = Q< Q< — i 
AjQ; = Q( + 1 — Q| 

A { Q ( = Q; + 1 - 2Q ( + Q, _ ] 

( v {A,) 2 Q ( - = Q, + 2 - 4 Q, + 1 + 6Qf - JQ, _ , + Q, _ t 
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Equivalent formulas arc used for differences in the n direction. 

A few details should be noted at this point. First, the sign in front of the artificial viscosity term be 
added to equation (7.5a) depends on the sign of the T term in the difference formula. For damping, t 
term must be negative when added to the right hand side of the equations (i.e., explicit artificial viscosi 
and positive when added to the left hand side (i.e., implicit artificial viscosity.) See Anderson, Tanncl 
and Pletcher (1984) for details. Second, the terms being added are differences only, and not finite differe 
approximations to derivatives. They arc therefore not divided by A£, etc. Third, the variables being 

A 

ferenced are Q, not Q. As noted by Pulliam (1986b), scaling the artificial viscosity terms by \/J makes th 
consistent with the form of the remaining terms in the equations. Fourth, the terms are also scaled by 
This makes the steady state solution independent of the time step size (Pulliam, 1986b). And finally, r 
that the fourth-order difference formula cannot be used at grid points adjacent to boundaries. At tl 
points, therefore, the appropriate fourth-order term in expression (8.1) is replaced by a second-order te 
Thus, for points adjacent to the £ = 0 and £ — 1 boundaries, — cSPAt^V^j^Q]/./ is replaced by 


+ 


J 


V ?A { Q 


( 


A similar expression is used at points adjacent to the rj = 0 and n = 1 boundaries. 

> t 

The implicit artificial viscosity is implemented by adding the following terms to the left hand side ol 
equations specified. 


to equation (7.5a) 

I 

to equation (7.5b) 

Note that the addition of the artificial viscosity terms, in effect, changes the original governing p: 
differential equations. At steady state, the difference equations with the artificial viscosity terms addec 
tually correspond to the following differential equations. 7 


c /At r a • .. 

— [VjAjfJAQ )] 

c /A t r a 

-~ 7 -[W 4 Q>] 


A A 

dE dF 
dt d*i 


ar„ eg 

d£ d n J 



gVQ) 


+ ( A >?) 2 


d\jQ) 

d n 2 




d\jQ) 

d? 


+ (An ) 4 


d\JQ) 

dn 


7 These equations represeni the use of the constant coefficient artificial viscosity model presented in this section, 
nonlinear coefficient model to be presented in Section 8.2 is more complicated, but the same principle appiic: 
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A 

The implicit terms do not appear, since they difference AQ, and in the steady form of the equations 
AQ = 0. The artificial viscosity terms do not represent anything physical. The coefficients should therefore 
be as small as possible, but still large enough to damp any instabilities. Although optimum values will vary 
from problem to problem, recommended levels arc = 0(1) and e, = 2c£ ) (Pulliam, 1986b). The recom- 
mended level for c£>, when used, is = 0(1). 

8.2 NONLINEAR COEFFICIENT ARTIFICIAL VISCOSITY 


The nonlinear coefficient artificial viscosity model is strictly explicit. Using the model as described by 
Pulliam (1986b), but in the current notation, the following terms are added to the right hand side of 
equation (7.5a). 


+ V, 


(T) l+ /(T)J(4%Q-4 4 V ! AQ).j 

(y) + (t)VV-<VAQ);: 

V ' j+\ V ' j_ 


(8.4) 


The difference operation A { V { A { Q is given by 

A « V {A{Q/ = Q / + 2 ~ 3 Qi + 1 + 3 Q/ - Q/ - 1 

In the expression (8.4), ^ is defined as 


^ = tx + ty 


(8.5) 


where tp x and i j/, are spectral radii defined by* 




+ v = - 


I U\ + ajt x 2 + { 
A£ 


V\+aJx 


2 

y 


nx + v y 2 


Ar, 


Here U and V are the contravariant velocities without metric normalization, defined by 


( 8 . 6 ) 


8 It should be noted that the grid increments Af and At? in these definitions do not appear in the corresponding 
formulas presented by Pulliam (1986b). This is because the grids used by Pulliam are constructed such that 
Ai = A/? = 1, while in Proteus A{ = l/fiVi — I) and Arj = 1/(^2 — 1). The definitions used here for i^ x and A, re- 
sult in an artificial viscosity level equivalent to that described by Pulliam. 
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U = Zt + Zx u + Zy v 

V = It + >!jc w + V 


and a = JyRT , the speed of sound. 

The parameters t m and c <4) are the second- and fourth-order artificial viscosity coefficients. Instcac 
being specified directly by the user, as they are in the constant coefficient model, in the nonlinear coeffic 
model they are a function of the pressure field. For the coefficients of the q direction differences, 

(eft = * 2 At max(a; + , , o t , <?,_,) (8 

(eft = max[0 , k 4 At - (eft] (8 

where 

— + 1 ~ ^ + Pi- l 

a ' _ A+i + 2 A + />,-! 

* i 

Similar formulas are used for the coefficients of the n direction differences. 

The parameter a is a pressure gradient scaling parameter that increases the amount of second-o 
smoothing relative to fourth-order smoothing near shock waves. The logic used to compute e w swit 
off the fourth-order smoothing when the second-order smoothing term is large. 

The parameters k 2 and k 4 are user-specified constants. Like the coefficients in the constant coeffic 
model, the optimum values will be problem-dependent, and are best chosen through experience. Cases ! 
been run with values of k 2 ranging from from 0.01 for flows without shocks to 0.1 for flows with she 
and 1C4 ranging from 0.0002 for flows computed with spatially constant second-order time differencin 
0.005 for flows computed with spatially varying first-order time differencing. Pulliam (1986b) ; 
kj = 0.25 and k 4 = 0.01 as typical values for an Euler analysis. 

Like the constant coefficient artificial viscosity model, the nonlinear coefficient model requires sp 
formulas near boundaries. To apply (8.4) at i= 2, is needed at i = 1. It is defined as 

(eft = k 2 At max(<7 2> cr,) 

With the above definition, applying (8.4) at /' = 2 and i = <V, — 1 requires a at /= 1 and i = AV The 
defined as 
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-Pt + 4pi - $Pi + 1p\ 
Pi + + 5Pi + 2 Pi 


°.v, 


~ P,V, - 3 + 4 Av, - 2 ~ SPN t - 1 + 2 P.V, 
Av, - 3 + 4 i°/V, - 2 + 5 Av, - 1 + 2 P/V, 


And, finally, applying (8.4) at i — 2 and i = /V, — 1 requires A^VjA^Q at i = 1 and i — — 1. 

numerous formulas that could be used. The ones currently in the Proteus code are 

= — Q 5 + 5Q 4 — 9Q 2 + 7Q 2 — 2Q ( 

A? V {A jQ^, _ , = _ 4 SQ^ _ 3 + 9Q iV( _ 2 - 7Q Nl - 1 + 2Q /V( 
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ABPetKi ix .Q 

Fluid Thermodynamic and Transport Property Equations 

The thermodynamic and transport equations used for the individual species and 
binary species pair combinations comprising a hydrogen- oxygen system is presented in 
this section. The equation types and specific coefficient values for the individual species 
thermodynamic and transport properties represent those used in the original input file of 
the CET89 computer code developed at the NASA Lewis Research Center (Reference 50) 
and are discussed in Section 0. 1 and 0.2, respectively. The binary species properties were 
based on recommendations made by Svehla (Reference 56) and are discussed in Section 
0.3. 


0.1. Thermodynamic Properties 

The fluid specific heat was represented by a piece-wise polynomial as described in 
Equation (0.1.1): 




( 0 . 1 . 1 ) 


— - — Si + a-i T + a-i + a A T^ + a« 


r. U 1 2 3 


l 4 a 5 

Likewise, the fluid enthalpy is described by the following polynomial: 

( 0 . 1 . 2 ) 


% = a, +^T + ^T 2 + ^T 3 + ^-T 4 + a 6 


RjT 


2 3 4 

and the fluid entropy is described in a similar manner: 


//~\ S i „ , „ t , 33 t 2 , a 4 t 3 , a 5 T 4 , „ 

(0.1.3) RT “ a l +a 2 ^ + + +a 7 

where the temperature in Equations (0.1.1) to (0.1.3) is in Kelvin. 

Low temperature and high temperature values of the coefficients presented in Equa- 
tions (0.1.1) to (0.1.3) were used, with the transition from the use of the low temperature 
to the high temperature coefficients occurring at a temperature of 1000 K. The coefficients 
used for all the species considered in the hydrogen-oxygen chemical rocket system are 
given in Tables 40 and 41. 


Table 40. Low temperature thermodynamic property coefficients. 



a t x 10 1 

a 2 x 10 2 

a 3 x 10 5 

a 4 x 10 8 

a 5 x 10 11 

a 6 x 10" 4 

a 7 x 10° 

H 

0.25000000 

0.0 

0.0 

0.0 

0.0 

2.5474390 

-0.45989841 

h 2 

0.29432327 

0.34815509 

-0.77713819 

0.74997496 

-0.25203379 

-0.097695413 

-1.8186137 

h 2 0 

0.41675564 

-0.18106868 

0.59450878 

-0.48670871 

0.15284144 

-3.0289546 

-0.73087997 

0 

0.30309401 

-0.22525853 

0.39824540 

-0.32604921 

0.10152035 

2.9136525 

2.6099341 

0 2 

0.37837136 

-0.30233634 

0.99492754 

-0.98189101 

0.33031826 

-0.10638107 

3.6416345 

OH 

0.38737299 

-0.13393773 

0.16348351 

-0.052133636 

0.0041826975 

0.35802349 

0.34202406 


Table 41. High temperature thermodynamic property coefficients. 



ajXlO -1 

a 2 x 10 3 

a 3 x 10 7 

a 4 x 10 10 

a 5 x 10 14 

a 6 x 10“ 4 

a ? x 10° 

H 

0.25000000 

0.0 

0.0 

0.0 

0.0 

2.5474390 

-0.45989841 

h 2 

0.30558123 

0.59740400 

-0.016747471 

-0.21247544 

0.25195487 

-0.086168476 

-1.7207073 

h 2 o 

0.26340654 

3.1121899 

-9.0278449 

1.2673054 

-0.69164732 

-2.9876258 

7.0823873 

O 

0.25342960 

-0.01247817 

-0.12562724 

0.06902986 

-0.063797098 

2.9231107 

4.9628592 

0 2 

0.36122139 

0.74853166 

-1.9820646 

0.33749007 

-0.23907374 

-0.11978151 

3.6703308 

OH 

0.28897814 

1.0005879 

-2.2048807 

0.20191288 

-0.03940983 

0.38857042 

5.5566427 
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0.2. Transport Properties 

The procedure for calculating the viscosity and thermal conductivity of a particular 
species is given in Equations (0.2.1) and (0.2.2): 


(0.2.1) [L = E T A exp (B/T + C/T 2 + D) 


( 0 . 2 . 2 ) 


k = ET A exp(B/T + C/T 2 + D) 


where the A, B, C, and D curve fit coefficients are unique for Equations (0.2.1) and 
(0.2.2) and have two sets of values that are valid in the two temperature intervals from 
300 to 1000 K and from 1000 to 5000 K. The values of these coefficients are given in 

_q 

Tables 42 to 45. In the viscosity equation the coefficient E is equal to 6.72 x 10 , 

while in the thermal conductivity equation E is equal to 4.01285 x 10 -4 . The units of 
the viscosity are lbm/(ft-sec), while the units of the thermal conductivity are ft-lbm/(sec 3 - 
R). The temperature to be used in Equations (0.2.1) and (0.2.2) is Kelvin. 

Table 42. Low temperature thermal conductivity coefficients. 

C x 10" 4 
-0.68759582 
-1.9701961 
0.55969886 
0.31668244 
0.2278508 
-0.37257802 


D 

4.3477961 

1.7545108 

-3.9259598 

1.78085307 

1.0050990 

-0.49894757 


B x 10 -2 


H 

0.58190587 

0.46941424 

h 2 

0.93724945 

1.9013311 

h 2 o 

1.5541443 

0.66106305 

o 

0.73824503 

0.11221345 

02 

0.81595343 

-0.34366856 

OH 

1.0657000 

0.45300526 
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Table 43. High temperature thermal conductivity coefficients. 



A 

B x 10~ 3 

C x 10 -5 

D 

H 

0.51631898 

-1.4613202 

7.1446141 

5.5877786 

H 2 

0.74368397 

-0.54941898 

2.5676376 

3.5553997 

H 2 0 

0.79349503 

-1.3340063 

3.7864327 

2.3591470 

0 

0.79819261 

0.17970493 

-0.52900889 

1.1797640 

°2 

0.80805788 

0.11982181 

-0.47335931 

0.9518913 

OH 

0.58415552 

-0.87533541 

2.0830503 

3.5371071 


Table 44. Low temperature viscosity coefficients. 



A 

BxlO -1 

C x 10~ 3 

D 

H 

0.58190587 

4.6941424 

-6.8759582 

0.91591909 

h 2 

0.68887644 

0.48727168 

-0.59565053 

0.55569577 

h 2 o 

0.7838778 

-38.260408 

49.040158 

0.85222785 

0 

0.73101989 

0.60468346 

3.5630372 

1.0955772 

°2 

0.61936357 

-4.4608607 

-1.3460714 

1.9597562 

OH 

0.78530133 

-16.524903 

12.621544 

0.69788972 


Table 45. High temperature viscosity coefficients. 



A 

B x 10 -1 

Cx 10 -4 

D 

H 

0.51631898 

-146.13202 

71.446141 

2.1559015 

H 2 

0.70504381 

3.6287686 

-0.72255550 

0.41921607 

h 2 o 

0.50714993 

-68.966913 

8.745475 

3.0285155 

O 

0.79832550 

18.039626 

-5.3243244 

0.51131026 

o 2 

0.63839563 

-0.12344438 

-2.2885810 

1.8056937 

OH 

0.58936635 

-36.223418 

2.3355306 

2.2363455 




































274 


0.3. Binary Diffusion Coefficients 

The binary diffusion coefficient model used in this study involves a polynomial curve 
fit of tabular predictions made by Svehla (Reference 61). The form of the curve fit model 
used in the Tethys model is shown in Equation (0.3.1): 

a, + a 2 T + ajT 2 + a 4 T 3 + a 5 T 4 + a 6 T 5 
(0.3.1) p y c“P 

where the units of the fluid static pressure, P, is atmospheres, the conversion constant C is 

—3 

equal to 1.076365 X 10 , and the temperature is expressed in Kelvin. The polynomial 
curve fit procedure was found to yield a sample coefficient of determination, r 2 (Reference 
60, page 405), of 1.000 for all of the tabular data, which indicates that the variation in the 
tabular data (which contained fourteen points ranging in temperature from 0 to 7000 K) is 
properly modelled by the polynomial curve to an accuracy of at least 0.05 %. 

The values of the coefficients used for every binary species pair considered in the 
Tethys model are presented in Table 46, with units of ft 2 /sec. 


Table 46. Binary diffusion coefficients. 


to.. 

* ii 



a 3 x 10 6 


■1 

EBBS 

h-h 2 

-6.5538 

30.157 


-4.5998 

5.7366 

-226.23 

h-h 2 o 

-6.6010 

-2.5389 

9.0676 

-6.8087 

3.1830 

1.2932 

H-0 

0.72345 

16.996 

6.9347 

-4.3702 

4.4986 

-183.79 

H-o 2 

-1.3452 

14.659 

7.9147 

-9.4989 

10.668 

-480.09 

H-OH 

-5.3002 

10.630 

8.9498 

-6.0576 

2.4253 

26.018 

h 2 -h 2 o 

-1.5355 

13.568 

7.0940 

-12.650 

18.090 

-947.24 

h 2 -o 

-3.6271 

21.194 

5.8545 

-2.4850 

1.9181 

-53.869 

h 2 -o 2 

-2.8965 

19.097 

5.6091 

-7.1199 

7.9971 

-377.46 

h 2 -oh 

-11.084 

28.537 

4.6334 

-1.9617 

1.4647 

-63.079 

h 2 o-o 

-74.065 

215.35 

27.764 

-68.597 

87.064 

-4202.8 

h 2 o-o 2 

-2.0372 

5.7462 

1.8950 

-2.3952 

2.6230 

-119.82 

h 2 o-oh 

-1.3893 

-2.8926 

2.8463 

-3.9455 

4.1437 

-177.03 

o-o 2 

-1.0001 

6.0491 

1.9629 

-1.8505 

1.9983 

-90.678 

O-OH 

-1.5942 

3.8541 

2.7872 

-3.0969 

3.2659 

-142.76 

o 2 -oh 

-2.0099 

4.9120 

1.9139 

-2.5363 

2.8288 

-129.95 
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